General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 







(NASA-CB-157329) ATMOSPHEBIC MODEL 
DEVELOPMENT IN SOPPOBT OF SEASAT. VOLUME 3 
PREDICTION MODELS Final Technical Report 
(Ocean Data Systems, Inc.) 165 p HC A08/MF 


Onclas 

25957 





OCEAIM DATA SVSTEIVIS, INC. 

2400 GARDEN ROAD. MONTEREY. CALIFORNIA 93940 o 408/649-1133 


Submitted to 

JET PROPULSION LABORATORY 
Pasadena, California 


ATMOSPHERIC MODEL DEVELOPMENT 
IN SUPPORT OP SEASAT 

VOLUME III “ PREDICTION MODELS 

Pinal Technical Report 


Prepared under 
JPL Contract Number 954668 
(Subcontract of NASA Contract Number NAS7-100) 


September 30, 1977 


Prepared by 
Robert E. Wellck 
OCEAN DATA SYSTEMS, INC. 
Monterey, California 


ABSTRACT 


As part of the SEASAT program of NASA/ a set of four 
hemispheric, atmospheric prediction models were developed 
for the Jet Propulsion Laboratory under Contract 954668 
CSubcontract under NASA Contract No. NAS7-100) . The 
descriptors applied to these four models, which use a polar 
stereographic grid in the horizontal and a sigma coordinate 
in the vertical , are : 

PECHCV - five sigma layers and a 63 x 63 horizontal 
grid, 

PECHFV - ten sigma layers and a 63 x 63 horizontal 
grid. 

PEPHCV - five sigma layers and a 187 x 187 horizontal 
grid. 

PEPHPV “ ten sigma layers and a 187 x 187 horizontal 
grid. 

Conservation forms of the difference equations based 
on the Arakawa technique are integrated using either a 
fifteen- or four-minute time step on a 381 km. or 127 km. 
grid (at 60°N) for the 63 x 63 or 187 x 187 models, respec- 
tively. Pressure-gradient-force terras are replaced by a 
single geopotential gradient on local pressure surfaces to 
reduce inconsistent truncation error (Kurihara modification) . 
Stress is applied at the lowest level. A nonlinear pressure 
smoother is used along with momentum and temperature diffusion 
to help control computational noise. The horizontal boundary 
conditions are insulated slippery walls. Centered time 

ii 


differencing with time averaging of the pressure-gradient- 
force term in the momentum equations is used. Robert time 
filtering of the temperature and moisture solutions is 
used for computationo.1 stability. 

The moisture and heat source/sink terms are modeled 
in a similar manner to those in the early Mints and Arakawa 
general circulation model. Terms representing evaporation 
and large-scale condensation, sensible heat exchange, para- 
meterized cumulus convection and precipitation, and solar 
and terrestrial radiation are included. Dry convective 
adjustment precludes hydrostatic instability. 

Initialization of the models is based on a pattern 
conservation technique to obtain objective analyses of the 
state parameter structure from the surface to 50 mb. 
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I. 


INTRODUCTION 


This document describes a set of four primitive- 
equation atmospheric forecast models for the Northern 
Hemisphere that were designed and developed under Contract 
954668 for the Jet Propulsion Laboratory, {Subcontract 
lander NASA Contract NAS7-100) as part of the SEASAT program 
of NASA. The models are hemispheric, using a polar stereo- 
graphic grid in the horizontal and a sigma coordinate in the 
vertical. The descriptors applied to the four models are: 

PECHCV - five sigma layers and a 63 x 63 
horizontal grid 

PECHFV - ten sigma layers and a 63 x 63 
horizontal grid 

PEFHCV - five sigma layers and a 187 x 187 
horizontal grid 

PEFHFV - ten sigma layers and a 187 x 187 
horizontal grid 

The 63 X 63 and 187 x 187 grids correspond to mesh lengths 
of 381 km and 127 km at 60 °N, respectively. 

There are two main sections describing these forecast 
models. First, Section II, the physical-mathematical 
description, gives the governing equations and describes 
the various physical/dynamical processes and their param- 
eterizations that are included. Second, Section III 
describes the program structure of the various models. 

In both sections, and particularly in Section II, the 


approach has been to describe the five-layer, 63 x 63 model 
as a baseline and then to describe the necessary modifications 
to either increase the vertical resolution to ten layers 
or to increase the horizontal resolution to 187 x 187, or 
both. 


II. PHYSICAL-MATHEMATICAL DESCRIPTION OP THE FORECAST 

MODELS 

This section gives the meteorological and mathematical 
description of the set of four forecast models; in other words, 
the governing equations , numerical techniques , algorithms 
and initialization procedure which constitute the theoretical 
basis of the models. 

The basic model equations, the primitive equations, 
are given in Section II-A in a form applicable to a polar 
stereographic map projection with an x,y coordinate system. 

Section II-E then describes the grid structure (hori- 
zontal and vertical) and boundary conditions for the hemis- 
pheric grid. Section II-C describes the heating and moisture 
source terms which include the solar and terrestrial radia- 
tion and sensible heating, large-scale condensation and 
convective adjustment, and the parameterized cumulus convec- 
tion and dry convective adjustment. Section II-D describes 
the numerical procedures used to approximate the solution of 
the primitive equations. 

Finally, Section II-E describes the methods used to 
create an initial state for a forecast from the data provided 
by the analysis models. 


A 


Basic Model Equations 


The governing partial differential equations for the 
forecast models, written in flux form, are similar to sets 
used by Smagorinsky ^ ai (1965) and Arakawa et ^ (1969) . 
These equations are listed below for an x-y map projection 
with map factor m and vertical coordinate, o. 

Momentum equation in the x direction: 
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3t 
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Momentum equation in the y direction: 
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Thermodynamic energy equation: 
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Moisture conservation equation; 


9 (irq) 
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Pressure tendency equation: 
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Hydrostatic equation: 

3j> _ RT 

3 a a 

In the above equations, 

u = X direction wind component 

V = y direction wind component 

TT = terrain pressure 

T = temperature 

q = vapor pressure 

w = “ - vertical velocity 

m = map factor 

f = Coriolis parameter 

=5 stress components 
D^,D^ - momentum diffusion terras 

Dy = temperature diffusion term 

H = diabatic heating term 

Q = moisture souroe/sink term 

$ = gz = geopotential 
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[II. 6] 
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B. 


Grid Structure and Boundary Conditions 


Integration of the forecast model equations is performed 
on Phillips (1957) sigma surfaces in which pressure, p, is 
normalized with the underlying terrain pressure, tt; that is: 



Figure II-l shows the sigma surfaces used in the forecast 
models. 

The horizontal wind components, u and v, temperature, 

T , and geopotentials , ® , are carried at the levels where 
cr = 0.9, 0.7, 0.5, 0.3 and 0.1 for the five-layer models, and 
at ff = 0.95, 0.85, 0.75, 0.65, 0.55, 0.45, 0.35, 0.25, 0.15 
and 0.05 for the ten-layer models. The moisture variable, 
g, is carried at the lower three and lower six of these sigma 
surfaces for the five- and ten-layer models, respectively. 

The vertical velocity, w, is calculated diagnostically from 
the continuity equation for the layer interfaces (a = 0.3, 

0.6, 0.4, and 0.2 for the five-layer models, and a = 0.9, 0.8, 
0.7, 0.6, 0.5, 0.4, 0.3, 0.2 and 0.1 for the ten-layer models) 
and is assumed to vanish at a = 1.0 and 0 = 0 . All variables 
are carried at each point (i,j) of the horizontal grid. 

The hemisphere is mapped onto a polar stereographic 
projection true at 60° north for which the map factor is 
given by: 


II -4 




1 + sin 60° 
^ 1 + sin (J) 


[II. 8] 


where $ is the north latitude. The square grid is centered 
about the North Pole (and extends approximately to the equa- 
tor, as shown in Figure II-2) . A 63 by 63 or 187 by 187 
horizontal grid system is used which gives a grid distance 
of Ax = Ay = 381 km and 127 km at 60 °N, respectively. 

The horizontal boundary conditions are similar to those 
used in the six-layer NMC forecast model described by 
Shtoman ^ ^ (1968) . The model is bounded by rigid, imper- 
meable vertical walls placed on the grid on the next to 
outermost row of grid points, as shown in Figure II-2. 

The boundary is both insulated and slippery; that is, no 
heat or momentum exchange is permitted across the wall, 
but a parallel flow is permitted along the wall. This means 
that along the lower boundary row where y is fixed, the 
boundary conditions can be written as : 

TTi = (Tru)^^ = 3 [H.9] 

Ti = T3; (irq)3^ = (Trq)3; $1 = O3 

where the subscripts refer to the first and third rows 
(row 2 being the boundary) . 



FIGURE II-2: NORTHERN HEMISPHERE POLAR STEREOGRAPHIC 

GRID (63 X 63). 


IS 
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Since the polar stereographic map projection becomes 
quite distorted in the corner areas of the grid and since 
these regions are not of meteorological significance for 
the forecast ^ the computed tendencies are truncated near 
the edge of the grid. From 20 °N to the Pole, full values 
of the computed tendencies are used. Below the tend- 

encies are set to zero. In the region between 5®N and 20°N 

the amount of the computed tendencies used varies according 

2 

to the relation [ (sin$ - sin5®)/{sin20‘’ - sin5°)3 . This 
procedure results in a persistence forecast in the outer- 
most areas of the grid and increases the stability of the 
computation near the edge of the grid. 
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C* Heating and Moisture Source Terms 


The diabetic heating rates in the forecast models are 
calculated in a manner similar to that of Mints and Arakawa, 
as outlined by Langlois and Kwok (1969) . The equation for 
the total heating can be written as 

® = ®SW **■ ^LW [II. 10] 


where Hg^ is the heating due to absorption of shortwave 
(solar) radiation; is the heating due to longwave (ter- 

restrial) radiation; Hj. is the contribution of sensible 
heat from the surface layer (cr=.9tocr=1.0 for the five- 
layer models and cr = 0.95 to c = 1.0 for the ten-layer models, 
a constant flux boundary layer) ; and is the release of 
latent heat resulting from phase changes of water substance, 
usually the latent heat of condensation. 

Section 1 describes the computation of the heating 
rates due to shortwave radiation, longwave radiation, 

and sensible heat, Hj.. Section 2 describes the large- 
scale condensation process which contributes to H^. Section 
3 describes the moist and dry convective adjustment proc- 
esses V7hich also contribute to 
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1. Solar and Terrestrial Radiation and Sensible Heating 


Diagnostic calculations of the heating? i.e., the tem- 
perature tendency, AT/At, represented by Hg^, and Hp, 

are computed every hour and saved for use over the next 
one-hour forecast period. The temperature tendencies are 
of the form 


AT _ g AF 

At C Ap 
P ^ 


[II.llJ 


where g is the acceleration of gravity? C is the specific 

IT 

A TT' 

heat of air at constant pressure and, ^ is the flux diver- 
gence. Flux computations depend on the vertical distribu- 
tion of temperature, moisture and pressure above each sur- 
face grid point? i.e., only vertical fluxes are computed. 

The computations of Hg^, and Hp are described in 

Subsections 2, 3 and 4, respectively. 


a. Solar Radiation 

The incoming solar radiation at the top of the model 
atmosphere is given by 


S = So COS y [11.12] 

2 

where Sg is the solar constant (taken to be 1.92 cal/cm min.) 
and y is the zenith angle which is a function of time of day. 


latitude and season. 

The zenith angle is given by 

cos y = sin^ sin6 + cos^ cosfi cosh [11.13] 

where <|i = latitude (0° at Equator, 90® at North Pole) 
h = hour angle 

6 = declination angle of sun 

and 

sin6 = 0.39785 sin [4 . 88578+0 . 0172*DAV 

+ 0.03342 sin (0.0172*DAY)-0. 001388 cos C0.0172*DAY) 

+ 0.000348 sin (0.0344*DAY) 

“ 0.000028 cos (0.0344*DAY) ] 

where DAY is the number of days measured from the first of 
the year. 

In order to compute the depletion of the solar beam, S 
is initially partitioned into two segments. The radiation 
with wavelength shorter than .9 microns (65.1% of the energy) 
is subject to scattering by the atmosphere, but not to absorp- 
tion. Wavelengths greater than .9 microns are subject only to 
absorption. The calculations that are outlined below follow 
the scheme constructed by Joseph (1966) . 

In order to compute the radiative fluxes, a simple form 
of middle cloud is assumed to exist between a = 0 . 4 and a- 0.8. 
An empirical formula for middle clouds; 

P -C 

CL = 1.3 (|-)^ - [0.46 + 0.6 ( \ ^ ) ] [1 - (|-)c3 [H.14] 

s 2 s 
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gives the amount of cloudiness# CL# in terms of the relative 
humidity and pressure at a = C. C is taken to be 0.7 for 
the five-layer models and 0,65 for the ten-layer models. 

= 400 and C 2 = 300 for the five- layer models# and = 300 
and C 2 = 350 for the ten-layer models. Limits of 1.0 and 0.0 
are set on the range of values given by equation [11.14]. 

The absorption of incoming solar radiation for the upper 
layer (a = 0 to a = 0.6) is given by; 

A 2 = 0.349S [ (l-CL)F^ + CL(l-a^)F2l [11.15] 

and for the lower layer (a = 0.6 to a = 1.0) by; 


= 0.349S [ (1-CL)F2 + CL(l-a^)F^] [11.16] 

where 
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|.pw(U) , 0. 303 
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[11.17] 

0.271 
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[11.18] 

0.271 
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0.271 
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[11.20] 


The symbol pw is the precipitable water given by: 


pw = i / 


g ^ Pi 


0. 622 


1000 




[ 11 . 21 ] 


and U refers to the amount of precipitable water above a = 0.6 


and T refers to the total amount of precipitable water. 

These absorption functions are from Manabe and Moller (1961) . 

The clouds’ albedo / a ^ is set at 0.5. Finally, the amount 

of solar radiation that reaches the surface, S , is needed 

y 

for the calculation of sensible heating over the land areas . 
Following Mints and Arakawa; 
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where a^, the ground surface albedo, is given as a function 
of cosy over the sea and as a function of snow cover over 
the land: a is the albedo of the cloudy sky and a is the 

SC 

albedo of the clear sky and is given by; 

P_ 


Ug = 0.085 - 0.10727 In [siny ( ] 


[11.23] 


b. Terrestrial Radiation 

The infrared fluxes atCT=1.0, 0.6 and 0 . 2 are computed 
in the manner suggested by Mints and Arakawa and are given by; 
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1.0 


* 4 4. 

^ 1.0 


[11.24] 


0.6 
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and 


*'0.2 ®’o .2 


0.8(1-CL)(F^ 

1 + I.75[pt«(L) + 1.2 pw(U)] 


57J16 


where St is the Stefan-Boltzmann constant and U and L refer to 
the amount of precipitable water above and below a - 0 . 6 / 
respectively. 

^1 0 ' ^0 6 ^0 2 fluxes computed by assuming 

the surface temperature is given by a downward extrapolated 

1 

temperature ^ 9 “ ^0 7 ^ five -layer models 

and 'fg = ^ ^^'^0 95 *" ^0 S5^ ten-layer models. The 

factor 0.8 is an empirical correction to account for the 
effect of CO^. 

In equation [11.24], T is given over the oceans by the 

9 

sea surface analysis, and over the land areas it is given 
by the solution of a heat balance equation explained in 

Section II-C-l-c on sensible heat transfer and evaporation. 

* 

To compute the approximate fluxes, F , the method devel- 
oped by Danard (1969) is used with an emissivity function 
for water vapor from Kuhn (1963). Thus, 

^1.0 “ ^ (1-CL) { -Tg+TQ^g$ (X1 )+Tq^^[^(X2)-<3>(X1)] } 

+ CL[T^_g-T^-TQ g)4>(X4)] ) [11.27] 

^0.6 " StCl-CL) [-Tgt(T^-T^^.g)^(Xl)+T^^^<f>{X3)3 [11.28] 
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and 

F* 2 = St C (1“CL) { -T^t(Tg-TQ_g) [${X2)“$(X3)] 

+ $(X3) > -CL ) Eli. 291 

where 


f(X) = O.OIOSX^ + 0.1675X + 0.542 

[11.30] 

XI = log^o 

[pw (L) ] 

[11.31] 

X2 = log^o 

Epw(L) + 1.2pw(U)] 

[11.32] 

X3 = log^Q 

[1.2pw(U)l 

[11.33] 

M 

X4 = 

Epw(l)] 

[11.34] 


where 1 denotes the lowest level . 

For these calculations, pw as given hy equation 


Eli. 21] is modified to: 


pw 
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j 0.622 I 


, p ,0.85 
^ 1000 ^ 


^273 j 0.5 


Eli. 35] 


which is the pressure- and temperature-corrected mass of 
water vapor in the atmospheric column. 

Finally, the fluxes are converted to heating rates by 


^LW2 


(F 


1.0 




, 3600 ^ 

\4*1000 ^ 


(°K/hr) 

‘ s^p 


Eli. 36] 


which is used below the a = 0.6 level. 


^LWl 



. 6 


F 


0.2 


) 


f 3600 , 
\ 4 . 1000 ^ 


^ (°K/hr) 

s p 


[11.37] 


which is used above the a = 0.6 level. 
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c. Sensible Heating and Evaporation in the 
Arakawa-Type Planetary Boundary Layer 


Sensible heating is computed as a function of the dif- 
ference of the surface temperature, T^., and the temperature 
of the air near the surface, T . Similarly, evaporation is 
computed as a function of the saturation specific humidity 
at the surface, Q , and the specific humidity near the sur- 
face, Q , except over land, where a Bowen ratio is computed 
to determine the flux of moisture. Difficulty arises in 
the forecast models since near-surface temperature, wind and. 
moisture are not predicted. To surmount this difficulty, the 
technique of Mints and Arakawa is used. 

The air near the surface; i.e., a = 1.0 to tr = 0.9, 
in the five -layer models and <r = 1.0 to a = 0-95 in the ten- 
layer models, is assumed to be a very thin boundary layer 
which does not absorb or store a significant amount of heat 
or moisture. The fluxes of sensible heat and moisture are, 
therefore, the same through the top and the bottom of this 
layer. Generally, for heat transfer, an equation of the 
form: 


r 


-K M 
H Az 


[11.38] 


is used where 0 is the potential temperature. 
III. 38] is modified by Arakawa to 

^ „ K r , '•”bL 

^ Pl.O ^p a*(0g^-0 ) “* Az 

fi + ^ 

^ Az •* 


Equation 
0 } 

-^] [11.39] 
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wheire AZ = - Zj^ q. In equation [11.39], the expression 

for Kg is allowed to vary with the stability of the boundary 

ic it 

layer. K , a and y empirical constants. The constant 

c 

Y allows for a somewhat more realistic variation of F between 
* c 

day and night. The subscript s refers to the value of a 
parameter extrapolated to the surface from the two lowest 
sigma levels. 

The sensible heat fl\ix through the bottom of the 
boundary layer is given by 

A 

»r = I’l.o 

where, for this computation, 

Vg = 0.8 [11.41] 

where Gu, the gustiness factor, is 2.2 m/sec. A drag 
coefficient = 0.002 is used. 

Over the oceans, T can easily be obtained from 
equations [11.39] and [11.40], since the sea surface temp- 
erature T , is specified. Knowing T , a heating rate due 
to sensible heating for the lowest sigma layer can be deter- 
mined. 

Over the land areas, T_ must be determined from a sur- 

g 

face energy balance equation: 

(l+r)Hj. - Sg + = 0 [11.42] 
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in which only the heat storage has been neglected. Data 
from Budyko (given in Sellers 1965) has been used to deter- 
mine r^ a Bowen ratio, as: 

r = 9.6 (sin<{))^ - 7.93 (sin(^) + 2.0 [11.43] 


in equation [11.42], the amount of solar radiation reaching 
the ground, S , is given by equation [11.22] , and the infra- 
red flux from the surface, p, is given by equation [11.24], 
Equation III. 42] is modified over ice-covered ocean to; 


tl+r)Hj, - Sg + = B(T^ - Tg) [11.44] 

2 

where B = 697.83 erg/cm sec “K is the ice conduction coef- 
ficient for three-meter -thick ice and T = 271. 2°K. 

w 

Evaporation over the ice-free oceans is determined in a 
similar manner to the heat flux. The equation resulting is: 


where 


E = 


Vs 

X . Vs ^ 
Az 




1 = 


K 


1 + 


* 

a (0 


BL“®s^ 


Az 


0.622 e 



[11.45] 


[11.46] 


[11.47] 


and 


= exp(Ae - 1^) [11.48] 

Ae = 21.656 and Be = 5418. 

The details of the derivation of equation [11.45] are given 
by Langlois and Kwok (1969) . 

Using the Bowen ratio from equation [11.43] allows the 
determination of the evaporation for land areas and ice- 
covered ocean as: 

E = H/r [11.49] 

The expression for r given by equation [11.43] has the effect 
of making the land surface at low latitudes a greater source 
of moisture than those at higher latitudes. 

2 . Large-Scale Condensation 

The large-scale condensation mechanism operates under 
the premise that synoptic scale motions do not contain super- 
saturated states. If a state of supersaturation is discovered 
(which may be due to excessive evaporation from the underlying 
surface or due to the effects of horizontal or vertical advec- 
tion of moisture ) , then sufficient moisture is removed (as 
precipitation with the latent heat of condensation) to achieve 
a saturated state at a warmer temperature. The technique 
suggested by Mints and Arakawa as outlined by Langlois and 
Kwok (1969) is followed in the forecast models. 
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A solution to the equation 


F(T) = Cp(T-T^) -* L[q^-qg(T)] = 0 [11.50] 

is sought; where q^ is the supersaturation mixing ratio 
associated with the temperature and is the saturation 
mixing ratio associated with the new temperature T, L is 
the latent heat of condensation which is taken to be 600 
cal/gm. 

Equation [11.50] is a non-linear equation. A second- 
order Newton-Raphson approximation to the solution is used, 
which is given by: 


T 


T 

o 


F(Tq) 


p"CT ) F(T) 

[1 + 5 -^ 

F' {T ) 


[11.51] 


Expressions for the derivatives of F(T^) must now be 
generated. F'(T) may be written as: 


F-CT) + 

P 


where 


9s ™ 



[11.521 


and e is given by equation [11.48]. Differentiating 

S 

q (T) and ignoring dp/dT gives: 

*" ‘■'o' =1 + 5- 7^^ 
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which may be further simplified to give 


P' (T^) 



[11.54] 


where : 


Y = 


BeqgtT,,) 

I 2 




The second derivative F' ' (T^) is given by: 


P' ' (T ) 
o 


§il 

Cp dT 


[11.55] 


[11.56] 


which, after considerable algebraic manipulation, can be 
written as: 

' ^^o^ " 4 I [11.57] 

p o o 


Using equation [11.51], the change in temperature may 
be obtained and from this, the change in mixing ratio can 
be computed by; 

C 

Aq = AT [11.58] 

The air is warmed by the amount AT and moisture is removed 
by the amount Aq from the supersaturated layer. 

Since moisture is forecast at the three lower levels 
in the five -layer models and at the lower six levels in the 
ten- layer models, a precipitation algorithm is required. 


This algorithm proceeds from the lowest level to the highest 
level/ and is outlined below. 

■ Examine the layer for supersaturation . If we are 
at the lowest sigma level and the layer is supersaturated/ 
then precipitation occurs; otherwise/ proceed to the next step. 
Examine the next lower layer for saturation. If 
the layer is saturated, then we allow the precipitation to 
pass through the layer to the surface where it is added to 
the total at that particular grid point; otherwise, proceed 
to the next step. 

• If the next lower layer was initially xinsaturated and 
the addition of the precipitation from the layer above results 
in supers aturation, we repeat the calculations summarized 

in equations [11.51] and [11.58]. The moisture removed is 
treated as precipitation; otherwise, proceed to the next step. 

* If the additional precipitation does not result in 
super saturation at the next lower level, we add the moisture 
to the layer and decrease the temperature accordingly. 

Finally, precipitation is given by the equation: 

Pr = inCur) Aq (cm) [11.59] 

where r = 1.0/0. 8 or 0.8/0. 6 for the levels in the five layer 
models, and r = 1.0/0. 9, 0.9/0. 8, 0.8/0. 7, 0.7/0. 6, or 0.S'/0-5 
for the levels in the ten layer models. 
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3 . Convective Adjustment 


In the atmosphere, the direct consequence of vertical 
static instability is the initiation of small-scale vertical 
convection which is very active in transferring heat and 
moisture from the surface through the lowest layers of 
the atmosphere. This process is very instrumental in main- 
taining the overall potential energy of the atmosphere against 
its destruction by radiative absorption and emission. 

In the forecast models, the primary purpose of the con-' 
vective adjustment mechanism is the removal of vertical in- 
stabilities, whether they be generated by a large-scale dy- 
namical process or by heat and moisture transfer from the 
surface to the lowest computational level of the model. If 
the vertical instabilities are allowed to grow, they will 
promote vertical convection whose scale is dictated only by 
the grid length in the model. This will lead to computational 
instability and destroy the solution. In order to control 
this type of instability and also to simulate vertical con- 
vection in the atmosphere, a convective adjustment process 
consisting of two parts, a parameterized cumulus convection 
and a dry convective adjustment, is performed at each time 
step. 


a. Parameterized Cumulus Convection 

Following a method devised by Arakav/a et (1969) , 
three types of parameterized cumulus convections are con- 
sidered. The first type is that occurring between levels 
2 and 3; the second type is that occurring between levels 
1 t 2 and 3 ; and the third type is that occurring between 
levels 1 and 2. Precipitation may result from types 1 and 
2f but not from type 3. Figure II-3 schematically illustrates 
these three types of parameterized convection. 

t 

Both the five-layer forecast models and the ten- 
layer forecast models use the same parameterized cumulus 
convection levels as shown in Figure II-3. In the case of 
the five -layer models, the convection levels coincide with 
sigma surfaces, while in the case of the ten-layer models, 
two adjacent sigma levels are combined to form a composite 
convection level with an averaged temperature and mixing 
ratio. Should cumulus convection then take place, the 
changes in moisture and temperature are apportioned between 
the two model layers on the basis of their mixing ratios. 

Three parameters used to indicate which type of 
convection can occur are the energy integrals 

S = Cp T + gz [11.60] 

h = Cp T + gz + Lg = S + Lq [11,61] 

E = C T + gz + Lq„ = S + Lq„ [11.62] 

]p So 
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TVP12 I 


TYPU II 


TYPE III 



FIGURE II-3: SCHEMATIC REPRESENTATION OF 

THREE TYPES OF PARAMETERIZED 
CUMULUS CONVECTION, WHERE n 
IS THE ENTRAINMENT COEFFICIENT 
AND C IS THE CONVECTIVE MASS 
FLUX. 


E is constant along a moist pseudo~adiabat and S is constant 
along a dry adiabat. Thus/ comparison of h and E for var- 
ious levels allows for a sufficient test of conditional 
instability. 

Type 1 convection will occur if 

> E 3 [11.63] 

If E 2 = , a moist adiabatic lapse rate would exist 

between levels 2 and 3. Since h 2 < E 2 / inequality [11.63] 
implies that a conditionally unstable lapse rate exists. 
Furthermore, inequality [11.63] can be rewritten as: 

^2 ^ Z ^^3 ■ [11.64] 

Type 2 convection will occur if 

> h 2 and h^^ > E^ [11.65] 

Conditional instability does not exist here between levels 
2 and 3, but does exist between levels 1 and 3. 

Finally, type 3 convection will occur if: 

E^ > hj^ > E 2 [11.66] 

Levels 1 and 2 are conditionally unstable; levels 2 and 3 
are not. 

The procedure now is to write down equations describing 
the budget for the moisture, q, and the dry static energy, S, 
for each type of convection. Before this can be done, the 


energy budget for the lov/er part of the cloud is considered, 
neglecting any accumulated storage of energy within the 
cloud . 

We use type 2 convection as an example. When the 

entrainment coefficient, rj / is greater than 1 , entrainment 

exists and the parameter E in the clouds, E , is given by 

c 

a mixing of h^ and h 2 in the lower part of the clouds, as 
follows : 

Ec = h 2 + I (h^ - h 2 ) [11.671 

because h is approximately conserved with respect to an 
air parcel. At level 3, 


whereas , for the environmental air : 


Therefore , 


where 


^3 = Vb = ^^3 + Lq^{T3) 


T 


C3 


- T. 


_1 

I+Y3 


(E -E_) 
c 3 


C 

P 


^3 



3 


[11.693 


[II . 703 


[11.713 


For type 1 and 3 convection, E^ is given by h 2 and h^, respec~ 
tively. From equation [11.703, one can obtain approximately: 


[11.72] 


Y3 <V^3> 


*^c3 I+Y 3 L 


For types 1 and 2, the convective terms in the moisture and 
dry static energy budget can be written in flux form as; 


aqj 

g 3t 

iPg sqj 


= nc [q„, + 


Y3 (E^-Ej) q2+q3 


s3 1+Yo 1* 


g 

at 

AP 2 

3q^ 

g 

at 

AP 3 

as^ 

g 

at 

APj 

as2 


_ r ^ 3“‘^2 , 1 ^^2 

~ ^ 2 n 2 ^ 


= C( 2 (type 2 only) 


E “E, S--S„ 

^ r c 3 , 3 2^ 

= nC It-tt: — + — o ] 


g at 


■I+Y3 2 


®3"®2 1 ®2'®1 

= nc r <-^^)l 


] [11.73] 


[11.74] 


2 n 2 


[11.75] 


[11.76] 


[11.77] 


and 


Ap, as, S_-S 

■'^’or type 3 convection: 

iP2 3q2 iPi Sqj 


g at g at 


^i-q. 


[11.78] 


[11.79] 


and 
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[11.80] 


&P2 _ Ap^ as^ _ S3^--S2 

g at “ g at ^ ^ 2 ^ 

where n is the entrainment parameter and C is the convective 
mass flux. 

Following Arakawa (1969) we now state that the 
solution of the above differential equations has a charac- 
teristic time scale, Xf such that: 


at 


(h-E) = - 


Ch-E) 

T 


[11.81] 


For type 1 convection, we can then obtain 


(hj-Ej) 


t>Q = - AP / 

gXi 2+Y3 I (h2'-E3) + (^^'^J) (S^-S.,)] 


3 


[11.82] 


where is given by the empirical relations; 


= 3600 [1+0.667 (h2"E3)] ^ 0*^ 


[11.83] 


and 


= 1200 I22(h2-E3) - 7] ^2“^3 - 


[11.84] 


For type 2 convection, a value of 1“K was used for 
therefore, we obtain: 


n = 


*' l "**2 


Ej+CpU+Yal-hj 


and 


Ap 

gx. 


(h,-E ) 
1 e 


^3-^2. . ,^1 ^2, 


[11.85] 


[11.86] 


[n (l*t*Y3) [c ^ 
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[11.87] 


where ^2 is given be the empirical relations 

T2 = 3600 [l+2.25(hj^-E3) ] ^1-23 <0,5 


and 


Tj = 3600 [15.75{hj^-E3)-5.75] > 0.5 [11.88] 


For type 3 convection, one can obtain: 


C i 

^3 


(h^-Ej) 




and 


[11.89] 


where T3 is given by the empirical relations: 

T3 = 3600 [1+0.667 {h^-E2>] ^i“E2 < 0.5 [11.90] 


T3 * 1200 [22(h^-E2) ~ 7] [11.91] 


Finally, precipitation from types 1 and 2 convection is 
given by 

Pr = ( l9M- - 60P .-.. 0 j^) p^ [11.92] 


b. Dry Convective Adjustment 

The dry convective adjustment process consists of 
checking the model atmosphere for static stability and mo- 
difying the temperatures when necessary to achieve static 


stability. This is accomplished by first computing the 
vertical profile of potential temperature, 0, by: 

th 

where it is the terrain pressure, it; the k sigma level 
and JC = = 0.2858. Starting at the bottom of the profile, 

we then test to see if 0, > 9, , +9^, where 0^ = 3®K. If 

this is tlie case, the layer is stable and no temperature 
modification is necessary. If it is not the case, then the ’ 
potential temperatures from layer 1 through k are modified 
to achieve a stable profile, using the average potential 
temperature and a fraction of the temperature difference 
0 per layer. (The fraction depends on the number of layers 

w 

whose temperatures are being modified. ) This process is 
continued until the entire profile has been checked for 
unconditional stability and modified, if necessary. The 
revised temperature profile is then recovered using equation 
[11.93], 
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D. NumGrical Techniques 


The numerical techniques used to approximate the solution 
of the forecast model equations are given in this section in 
two parts. First, the spatial finite differencing tech- 
niques are described. Then a description of the time dif- 
ferencing technique is given which includes an explanation 
of the pressure-gradient-force averaging technique, Robert 
time filtering and the non-linear pressure smoother. 


1. Spacial Finite Differencing 


Energy conserving (in space) difference equations 
based on the Arakawa technique (1966) are used to approxi- 
mate the spacial derivatives. Thus, the finite difference 
equation for the momentum equation in the x direction is 
given by: 


3 (ttu) 
9t 


n 

[-L(u)* + f(ir*v*) 


, n r.n-1 


t D 


n-1 

u 


[11.94] 


and in the y direction by: 


[-L(v)» - fCx.u,) - [11.951 

where the subscript * denotes the point (i,j,k); the super- 
scripts denote the time level; L is the Arakawa advective 


9 (ttv) 
9t 
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si 


t 


operator; f is the Coriolis parameter; the A notation 
refers to the Kurihara pressure gradient; and the bar 
denotes time averaging as explained in the next section. 
F and F are the frictional terms; D and D are the 

X y U, V 

lateral diffusion terms. 

The Arakawa advective operator is given by : 


ilL • • 




III. 961 


TT . . 


t«k - "k-i'Wi'hj 


where a = ~ 
m 



m = map factor 
d = grid mesh size 


The pressure gradient terms ^ and ^ are computed 
using a type of modification suggested by Kurihara (1968) . 
Instead of evaluating both the gradient of the geopotential 
distribution on a sigma surface and the terrain pressure 
gradient (invariant with height) as would be required by 
equations [II. 1] and [II. 2], the value of the geopotential 
is interpolated at neighboring grid point locations to the 
same pressure surface as at the computation pointy allowing 
computation of the geopotential gradients on pressure surfaces 


The interpolation formula used is a form of: 


9 $ 

aino 


-RT 


tll‘973 


This procedure of differencing $ on local pressure surfaces 
reduces the inconsistent truncation error over high, irreg- 
ular terrain. 

Frictional dissipation is accounted for in the gross 
boundary layer (a = 1.0 - 0.8 for the five-layer models and 
o = 1.0 - 0.9 for the ten-layer models) only. The stress is' 
assumed to vanish at the top of this layer. The stress 
terms are the same as those given by Shuman and Hovermale 
C1968) , except that the surface wind speed is obtained by 
extrapolation from the boundary layer sigma surface. Thus, 
the friction in the x direction is given by; 


,n-l 


AuR ''s '^T 


[11.98] 


and in the y direction by 


,n-l 


r.n-1 /TTVn 

AaR '^s '^T ^ 


n-1 


BL 


[11.99] 


where Vg = 0.3[Ugj^ + 


[II. 100] 


and Cj^ = 0.0015 for sea level and 0.0025 for non-sea level 
terrain. 
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The lateral diffusion terms are given by 


m? 

D, = D (irV^u)5~^ 

d 


D = D (irV^v)5"^ 

V 


[II. 101] 


[ 11 . 102 ] 


where the laplacian is 




[11.103] 


and the diffusion coefficient, D, varies with latitude according 

2 

to the relation K (1-sin cj)) with 0.9 taken as a maximum for 

g 

sin^. K is taken to be 10 for the 63 x 63 models and 
4 X 10^ for the 187 x 187 models. 

The change in the terrain pressure is computed from 


an approximate form of the continuity equation as: 


[11.104] 


r, K ^2 

3TT ^ 1 V* ^i 3 

9t ^ "K ~2t‘ ^“i+l,j,k"°‘i-l,j,k'^^i,j+l,k"^i,j-l,k^ 
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The thermodynamic energy equation is approximated by: 


9 (ITT) 
9t 

n 

Oj 

* 

- L(T)J 

[11.1051 

cp ' 

( 

2CTj^ ^^k‘**'^k-l^ i, j ^ ^&t^i,j 
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All 


"23“ ’^^i, j,k^'^i+l,j,k ^i-l,j,k^ 


^i, j ,k^'^i, j+l,k“'^i, j-a,k^ ^ ^ 


n 


, „ ^ n— 1 , „n“l 
+ H^'FT . . + D„ 

* 1,: T 


where is given by an expression similar to those for 
and (Equations II-lOl and 11-102) . 

Vertical motions for non-material surfaces are obtained 
from a form of the continuity equation; 


'^i,j,k+l = 


Aa ( ,3ir^ 

7T. . ^ '3t^i,j 


[11.106] 


^ 2d^ ^^*^1+1 “i-1^ j^k'^^^j+l ^j-l^i,k^ ^ 


The hydrostatic equation is given at the lower level by: 




[11.107] 
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where is extrapolated from the lowest two levels by 
1.25Tj^ - O. 25 T 2 and is the terrain geopotential. For 
subsequent levels 


k+1 


= »k - 1 


+ ’’k+i' 


[11.108] 


Finally, the moisture equation is approximated by 


3 (TTq) 
3t 



“ L(q) 


n 

* 






n-1 

i. j 


[11.109] 


2. Temporal Finite Differencing 

A centered, or leap-frog, difference scheme is used 
to approximate the time derivative in the forecast model 
equations ; that is 

2^t [II. 110] 

In order to lengthen the permissible time step used in 
the computation, a numerical technique known as pressure- 
gradient-force time averaging is used in conjunction with 
the centered time differencing. This procedure was suggested 
by Shuman (1971) and has been successfully used at the 
National Meteorological Center. Pressure-gradient-force 
time averaging requires arranging the integration of the 
primitive equations into two parts. First, the pressure 
tendency, thermodynamic energy and moisture equations, in- 
cluding the diabatic terms, are integrated forward to the 


3A 

3t 
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(n+1) time level. Then the momentiom equations are integrated 
forward in time to the (n+1) time levels but the pressure 
gradient force terms (refer to Equations [11.941 and [11.95]) 
are time-averaged according to the relation; 

= (l-a)PGF^ + + PCF^"^^) [II. Ill] 

For stability, 0<a<0.5, and it can be shown that when a=0.5, 
a time step twice as long as when a— 0 can be used. However, 
a weak computational instability is present and a value of a 
slightly less than one-half should be used. The value 
a=0.49 is used in the forecast models. 

A centered difference scheme is neutrally stable in 
that there is no damping, only phase shift, of the various 
frequency components present in the solution of the difference 
equations. In order to achieve damping of the shorter wave 
length components to avoid non-linear instability and, in 
particular, to damp the computational mode, the simple leap- 
frog scheme is modified slightly to time -filter the solutions 
of the thermodynamic and moisture equations . This procedure 
is known as time filtering and was originally used by Robert 
(1966) . 

Assume that the partial differential equation that we 
wish to approximate is 

|^=F(A) [11.112] 


11-38 


Then, using the leap-frog difference scheme, we have 

+ 2At F(A) [11.113] 

This step is modified by Robert to become 

+ 2At P(A,A*) [11.114] 

A^ = (1-u) (A’)^^ + f [A^"^ + (A*)”'^^] 0<a<l [11.115] 

where we time-filter the solution. It can be shown that this 
procedure adds damping to the solution and heavily damps 
the computational mode that is present in the solution of 
the finite difference equations. 

If ojAt<l is the linear stability criteria for the 
simple leap-frog difference scheme, then the addition of 
Robert time filtering will change the stability criteria to: 

uAt < [11.116] 

A value of a=0.1 is used to time-filter the solution of the 
thermodynamic and moisture equations. No time filtering of 
the solutions of the momentum equations is done since lateral 
diffusion terms are included in these equations. 

The terrain pressure tendency equation [11.104] is 
essentially an ordinary differential equation? thus, rather 
than using Robert time filtering of the solution to control 
small-scale noise and solution separation, a non-linear 
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smoothing operator is applied every time step. This smoothing 
operator has been described by Oliger and Wellck (1970) and 
is given by: 


APS = 


6Ps 




' max 


i+1/: 


1 „n+l _n I /t,n+l T^n ^ 

- 1PS» - P=i-I,jl (PS. - P®i-l,j> 


^ |p§i,3^i - pSt"i (Pii,3^i - psr> 


n+l, 


Ps 


n+l «n 




n+l 




[11.117] 


where ISPsl ^ is the maximtim absolute difference between 
Ps at the point (i/j) and Ps at any of the four immediately 
neighboring points, and s is the smoothing coefficient which 
has been determined to be of the order of 0.01 for control 
of computational noise. 

The smoothed surface pressure is then given by: 

Ps (smoothed) = Ps + APs [11.118] 


The terrain pressure is reduced to sea level pressure and 
vice versa through the relations 

Ps = ir^ + Att^ [11.119] 


and 


Att^ = Air, 


[1+(J- - 1) 
^0 


,tt^-700, 

^ 300 ^ ^ 


[ 11 . 120 ] 


where tTq is the initial terrain pressure and Attq is the 
initial difference between sea level and terrain pressure. 

These relations are also used to create surface pressure 
fields for output. 

Due to the interaction of these various computational 
devices described above (pressure-gradient-force time av- 
eraging,. Robert time filtering, lateral diffusion and non- 
linear pressure smoothing) a doubling of the permissible 
time step is not possible. However, a fifteen-minute time 
step for the 63 by 63 models as opposed to a ten-minute 
time step, and a four-minute time step for the 187 by 187 
models is possible which results in a considerable saving 
of computer time over simple centered time differencing. 

Finally, the integrations are restarted after each 
output cycle (normally every twelve model- hours) using an 
Euler-backward or Mats\mo time step in the 63 x 63 models. 
Provision is made to perform an Euler-backward time step 
(which is highly dissipative) at specified intervals in 
tlie 187 X 187 models. This integration step is given by: 

(A')^'*'^ = + At P[A^3 [11.121] 

and 


« A^ + At F[(A’)^’^^] 


[ 11 . 122 ] 


E, 


Initialization 


In order to be able to forecast the meteorological 
conditions at some time in the future using a primitive 
equation model the conditions at the present time must 
be knovm; that is, an initial state must be specified for 
the variables forecast by the PEM* This initialization 
process is accomplished by using analyzed data produced by 
a set of analysis programs and synthesizing this data into 
a complete specification of the initial state for the PEM. 

The following sections describe the method of inter- 
polating the analyzed data fields to the forecast model 
sigma surfaces. The specification of various other required 
data fields is also described. 

1. Constant Fields 

Constant data fields are of two types — those that 
are truly constants (sine of latitude, map factor, tendency 
truncation coefficient and terrain height) and those that, 
while not constants , are held constant during the course of 
a forecast (sea surface temperature, land-sea- ice indicator 
and albedo) . 

For a polar stereographic grid, true at 60 and centered 
about the North Pole, the sine of latitude is given by: 

= STP ^It 111 . 123 } 
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where 


I 

1 

1 


STP = (31.205)^ for the 63 by 63 grid 
STP = (93.615)^ for the 187 by 187 grid 

and 

ST = (I-POL) ^ + ( J-POL) ^ 

POL = 32 for the 63 by 63 grid 
POL = 94 for the 187 by 187 grid 


The map factor is given by 


1 + sin 60° 
1 + sin (p 


[11.124] 


The tendency truncation coefficient used to reduce the 
tendencxes of the variables in the outer regions o'*^ the 
grid is given by: 


re = 1.0 for ^ > 20° 


[11.125] 


= 0.0 for (|). < 5° 

and varies as [(sin^ - sin5°) / (sin20° - sin5°)]^ for 4) between 
5° and 20°. 

The terrain height is derived from Scripps Institute of 
Oceanography data and has been are a- averaged and gradient- 
limited, That is, the terrain gradient at any grid point is 
constrained to be less than 2000 m per 381 km, the basic 
63 by 63 grid size. In addition, the terrain is multiplied 
by the tendency truncator thereby tapering it towards sea 
level around the edges of the grid. 
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The land- sea- ice indicator is also derived from the 
geography and varies seasonally as ice masses expaa^d and 
contract . 

The albedo, ALB, of the earth's surface is obtained 
from the seasonal albedo charts of the northern hemisphere 
developed by Posey and Clapp. The open sea albedo is 
modified in the heating computation as a function of the 
zenith angle of the sun according to the relation; 

ALB' = ALB +0.54 (0.7-cosy) fll.126] 

The sea surface temperature is one of the analyzed data 
fields produced by the pattei-n conservation technique analysis 
programs . 


2. Terrain Pressure 

The terrain pressure, tr, is determined from the ana- 
lyzed sea level pressure, P^, the terrain height, and the 

heights of the two analyzed pressure surfaces bounding the 
terrain height. That is, if the terrain height lies be- 
tween and Zy, which correspond to the heights of the pres- 
sure surfaces, P^^ and P^^, then the terrain pressure is given 
by 

Pl ^ 2 -Z 

IT = P . (^) where X = ■■ [11.127] 

U Py Zy-Zj^ 

The difference betvi/een the sea level pressure, P_, and 
the terrain pressure, fr, is also computed and stored for use 
in reducing the forecast terrain pressure to give a forecast 
sea level pressure. 


11-44 


3. Temperature 


The sigma surface temperatures are determined by loga- 
rithmic interpolation between the analyzed pressure surfaces. 
The terrain pressure, it, and the temperatures of the analyzed 
pressure surfaces, are used to determine the temperatures of 
the sigma surface from the following relation: 


InCP^) “ In (air) 

^ •^ln(Pj^) - In (P^) ^ 


til. 128] 


That is, the pressures, and P^, are found which bound 
the sigma surface, cttt, and along with the corresponding 
temperatures, T^ and Tt,, are used in equation [11.128]. 

JLj u 

4- Geopotentials 

Having computed a set of sigma surface temperatures , and 
having obtained the terrain geopotential, the hydrostatic 
equation can be used to compute an initial set of geopoten- 
tials. That is, equations [11.107] and [11.108] are used to 
compute the sigma surface geopotentials. 

5. Winds 

The initial winds are obtained from an analysis pro- 
cedure which gives analyzed winds on pressure surfaces from 
1000 to 100 mb, and are computed geostrophically on the 
50 mb pressure surface. The 63 x 63 models use wind data 
on eleven pressure surfaces: 950, 850, 700, 500, 400, 300, 
250, 200, 150, 100 and 50 mb. The 187 x 187 models use 
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wind data on all thirteen of the pressure surfaces. The 
initial winds on pressure surfaces are then vertically 
interpolated linearly in the logarithm of pressure to the 
model sigma surfaces using equation [11.128]. 

6. Moisture 

An important part of the initialization process is the 
specification of the initial moisture distribution. The 
forecast models carry moisture at the lower three levels in 
the five-layer models and at the lower six levels in the ten- 
layer models and would, thus, require a large amount of ana- 
lyzed moisture data for initialization. Since moisture data 
are virtually non-existent over the oceans, a parameterized 
approach to the specification of the initial moisture dis- 
tribution has been taken. 

Studies by Kesel and Lewit (1974) showed that existing 
analysis/procedures for moisture initialization do not 
necessarily result in an optimum moisture distribution. It 
was further shown that improved forecasts could result from 
moisture initialization based upon the relative vorticity 
distribution at each model level. 

ODSI’s studies have shown that the initialization of 
moisture has a significant effect on deepening lows, but 
has little effect on filling lows or highs. It was shown 
that by basing the initial relative humidity on the intensity 
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of the geostrophic relative vorticity, specific regions of 
properly posed relative humidity values will more likely co- 
incide with actively developing systems . 

Therefore, a procedure based on the intensity of the 
geostrophic relative vorticity is used to define the initial 
moisture distributions on the 1,000, 900, 700, 500 and 300 mb 
pressure surfaces. These values are then interpolated linearly 
in the logarithm of pressure using equation [11.128] to the 
sigma surfaces. This procedure for each of the five pressure 
surfaces can be outlined as follows : 


a. Compute the geostrophic relative vorticity 
for each point in the field using the relation 



[11-129] 


where Z is the height in meters of the pressure surface and 
f is the Coriolis parameter. 

b. Compute the mean and standard deviation of 
the relative vorticity for the field using the relations 


M - 


1 

N 


N 

2 



[11.130] 


SD = [i 1 j 


[11.131] 


c. Normalize the relative vorticity field according 
to the relation 


RV . . -M 


[11.132] 
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d. Compute the relative humidity for each point 
based on the normalized relative vorticity values using the 
empirical relation 




[11.1331 


with the upper limit of relative htamidity set at 100% and 
the lower limit set at 20%. The parameter REL depends on 
the pressure surface and is given by 0.7, 0.6, 0.5, 0.4 
and 0.3 for 1000 through 300 mb, respectively. 

e. Modify the relative htimidities computed using 
equation [11.133] according to whether they are over ocean 
or land or are in thesntial lows (summers) . If the point is 
over the ocean, then the computed humidity is unchanged, 
while if the point is over land the computed humidity is 
multiplied by 0.8. Further, if the point is in a thermal 
low, then the computed humidity is multiplied by 0.5. A 
thermal low is defined as existing when the day as measured 
from January 1 is between 150 and 250 (summer) , and either 
the 1000 mb - 500 mb thickness is greater than 5790 meters 
or the 500 mb temperature is greater than -5°C. 
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III. COMPUTER PROGRAM DESCRIPTION OF THE FORECAST MODELS 

The designer of any large, time-consinning computer 
code is faced with n\imerous problems in trying to design an 
efficient and readable program. First, there are the prob- 
lem constraints on how many grid points are desired, how 
many variables are to be carried at each point, exactly 
what must be computed, what are the inputs and what are the 
desired outputs. Second, there are the hardware constraints 
— the size of main memory, the type of se :>ndary storage 
available (rotating or extended memory) , the number and 
characteristics of the I/O channels and the speed of the CPU 
Third, there are the software system constraints — exactly 
what will the system allow one to do, what sort of high- ! 
level language is available and what are its characteristics 
Finally, there is the requirement that the resulting code 
should be efficient, readable and easily modified, both from 
the standpoint of computational changes and portability 
between computers . 

What follows is a description of the realization of 
these various design objectives, including the various com- 
promises that were required. There are four forecast models 
described that range in data base size from approximately 
200,000 to 3,000,000 words and span roughly a factor of 
fifty in computational requirements. The forecast models 
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use two basic data flow structures and are grouped as 
follows: the five- and ten-layer/ coarse horizontal grid 

models r PECHCV and PECHPV, and the five- and ten-layer 
fine horizontal grid models, PEFHCV and PEPHPV. 

During the course of developing these forecast models, 
three different computer systems were utilized: a GDC 6500, 

a CDC CYBER 175 and a CDC CYBER 76, all using CDC software. 
Thus , an attempt at portability was required and largely 
realized, although the forecast models are in no sense 
totally portable between all computer systems with FOR- 
TRAN IV as their high-level language. 

It is not the intent of this section to give complete 
detail down to the definition of the last local variable in 
each subprogram, but to describe the data structure of each 
program and to describe the code structure with sufficient 
detail to be easily understandable without overwhelming the 
reader with trivia. Also, since there is a large amount of 
duplication of code between the four forecast models, dif- 
ferentiation between programs in the explanations which 
follow is only made when necessary. 

Section A explains the data structures of the four 
forecast models. Section B describes the program structures 
of models PECHCV and PECHFV and Section C describes the 
program structures of models PEFHCV and PEPHPV. Included 
in these sections are block diagrams, definitions of the 
various common blocks and functional descriptions of each of 
the subprograms. 
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A. 


Data Structure 


One of the most critical factors involved in attempting 
to design an efficient computer program is an efficient 
data storage scheme, since a primitive equation forecast 
model (PEM) requires the handling of large amounts of data. 
Section 1 gives a general picture of the data storage 
philosophy used in the models. Section 2 then discusses 
the details of the two versions of data storage structure 
used in the models giving diagrams, tables and illustrative 
examples of the working of this method of data storage. 

1. Data Storage - Philosophy 

Most PE models numerically approximate the solution 
of a complicated set of partial differential equations — 
the primitive equations — through the use of finite differ- 
ences. While these PEM are three-dimensional in space, the 
treatment of the third space dimension, the vertical dimen- 
sion, is quite different from the two horizontal dimensions. 
The use of the hydrostatic assumption to render the numerical 
solution feasible and vertical colximn computations in the 
heating packages immediately offer themselves as obvious 
examples of the difference between the horizontal and the 
vertical dimensions. 

If one assumes some sort of grid system, for 
example, either a polar stereographic projection or a 
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spherical polar grid, one then has an array of data, say MxN 
in the horizontal with K vertical levels. Since there will 
be at least five variables at every point, the total number 
of pieces of information can be quite large, on the order 
of millions of words. Except for the very largest computer 
memories, this vast amount of data, together with the nec- 
essary computer code, cannot be contained in memory. Thus, 
one is faced with the problem of devising some sort of 
efficient storage scheme for the data. 

Meteorologists are accustomed to looking at weather 
data on pseudo-horizontal surfaces; for example, the geopo- 
tential heights at 500 mb. Very few outputs from a PEM are 
presented as vertical slices. This, of course, is physically 
quite reasonable, since the scales of motion are quite 
different between the horizontal and the vertical, as are 
the grid distances in a PEM. However, horizontal planes 
are not the best way to organize the data storage for com- 
putation in a PEM. 

Since it has been assumed that the numerical method 
used to approximate the solution of the primitive equations 
is finite differencing, it is easily seen that an entire 
horizontal plane is not necessarily needed in memory simul- 
taneously. The second order finite difference approximations 
for horizontal derivatives used for the forecast models 
are given by 
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3x 


ir j/k 


[III.l] 


2 Ax "" 


and 


9u 

9y 


'if jrk 


2Ay “ “i,j-l^k 


[III. 2] 


If one imagines the j subscript being held fixed and i 
varying, then there are only three j rows involved in 
the computation along a single j row, the one being com- 
puted, the one below and the one above. In addition, since 
the heating package computations are all performed on ver- 
tical columns, it would seem reasonable to use a data 
storage scheme such as the one shown in Figure III-l. 

The scheme shown in Figure III-l has the computa- 
tional data broken down into vertical slices which are MxK. 
There are N of these slices. It has been seen that only 
three of these slices are needed at any one time to compute 
along a given j row. Thus, if one were to construct a file 
with N records on it, each containing one vertical slice, 
one could perform a forward integration of the primitive 
equations by sweeping this file in sequence. Of course, a 
new file should be constructed as this file is swept, in 
order to have data for the next step. The details of such 
a file scheme will be explained in the next section. 

A PEM consists of three basic sections, the ini- 
tialization section, the computational section which 
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FIGUP3 


III-l: ROW OR VERTICAL SLICE 

DATA STORAGE SCHEME 
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actually integrates the primitive equations, and an out- 
put section. By far the most computer time is spent in 
the computational section; thus, it seems reasonable to 
arrange the data storage for maximum convenience of this 
section of the code. This hypothetical PEM would take, 
as input, data on horizontal planes (pressure surfac.?), 
since this is the way that analyses are performed. It 
would then interpolate this data to the sigma surfaces and 
finally sort out this data into vertical slices for use by 
the computational section. The computational section would 
then integrate the data forward in time to obtain a fore- 
cast and pass it to the output section. The output section 
will then sort out the data into horizontal planes and 
interpolate, smooth, etc., for viewing by the users of the 
forecast. 

An obvious question that arises is , does this 
method of data storage result in an efficient PEM from the 
standpoint of central memory and computational speed? The 
answer is that it does and that this method of data storage 
is in use at numerous meteorological centers. As will be 
explained in the next section, rotating storage can be used 
quite effectively and more efficient utilization of central 
memory can be obtained when using the vertical slice storage 
technique. 


2 . 


Data Storage -> Realization 


The concept of a data file consisting of N records, 
each representing a vertical slice of the three-dimensional 
grid system, was introduced in the previous section. If 
the grid is MxNxK, then there will be N records on the file. 
Each record will hold data for an MxK vertical slice. It 
was also pointed out that, in order to perform computations 
on the j row, only the j-1 and j+1 rows are needed in addi- 
tion to the j row. In this section, the file structure, 
memory storage requirements and data flow will be examined 
in detail. 

For the sake of clarity, let us call the existing 
data file lOLD and the data file created by one sweep of 
the data base INBW. If it is assumed that centered time 
differencing is used, i.e.. 


3A 

3t 
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n+1 
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and that the numerical technique of pressure-gradient— force 
averaging is used, then file lOLD will contain information 
at the n-1 and n time levels. After the first sweep of lOLD, 
INEW will contain partially updat">d information at the n and 
n+1 time levels (temperatures, moistures and geopotentials) 
along with the non- updated information. At the end of the 
sweep, the file INEW becomes the file lOLD and a second 
sweep of the data is made to integrate the momentum equations 
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thereby completing the time step. The file INEW now con- 
tains all updated data at time levels n and n+1. This 
process is repeated until the desired forecast time is 
reached. 

In order to achieve more efficient use of central 
memory and to reduce the size of the data files f variables 
with time levels are "time packed" on a 2:1 basis. That 
is, the n time level occupies the upper thirty bits of the 
word and the n-1 time level occupies the lower thirty bits 
of the word. Certain other infrequently referenced vari- 
ables are also packed on a 2:1 basis. A set of constructs 
which take the form of FORTRAN statement functions have 
been designed to be able to manipulate these packed data 
words. These statement functions are listed below with a 
brief explanation of their function; 

retrieving the n-1 time level 
or lower value - 

SNM{V) = SHIFT (V, 30) [III. 4] 

pack two values into a word - 

PACK(VNP,VN) = (VNP.AND.77. . ,70. . .OB) .OR. 

(SHIFT(VN,-30) .AND.O. . .07. . .7B) [111.5] 

store a value into the top half 
of a packed word - 


STOR(VNP,VN) = (VNP.AND.77. . .70. . .OB) .OR. 

(VN.AND.O. . .07. . .7B) [III. 6] 


since this "time packing" method relies on the CDC 
FTN compiler intrinsic function SHIFT, the forecast models 
are not portable to other compilers. However, this was thought 
to be worthwhile since it affords a considerable compression 
in data file size and has a mini.mal effect on the computation 
time (approximately a 3% increase) . 
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a. 


Data Storage - PECHCV and PECHFV 


The concept of an old data file being read and a 
new data file being created on each sweep of the data bape 
was presented in the previous section. In addition, the 
idea of "time packing" and the necessary methods of acces- 
sing the packed data were also presented. The next item 
to be considered is an efficient means of flowing the data 
to and from central memory in order to implement the com- 
putation, 

A diagram of the necessary central memory storage 
buffers for this data storage scheme to work is given in 
Figure II I- 2. Two buffer areas are required, AO for the 
three working rows and the row being input, and AN for the 
row being created and the row being output. It has been 
assumed that the computer hardware will allow simultaneous 
I/O and computation- buffered computation; thus, the row JIN 
in array AO and the row JM in array AN. Further, since it 
is time consuming to move data about in memory, the buffers 
AO and AN are circular buffers, with the pointers redefined 
as the computation proceeds. 

One more construct is needed before a given piece 
of data can be accessed in the arrays AO and AN. This is 
a method of pointing to the particular variable that is 
required on the particular row that is desired. For com- 
putational convenience, two subscript functions have been 


Input Row, JIN 


from lOLD 


JP Row 


J Row {n-l,n levels) 


JM Row 


AO Array- 
Pointers JMN 
JN, JPN, and 
JIN 


J Row {n,n+l levels) 


Output Row r JM 


to INEW 


An Array 

Pointers JMNP and JNP 


FIGURE III-2; CENTRAL MEMORY STORAGE 

BUFFERS REQUIRED FOR ROW 
STORAGE SCHEME USED IN 
MODELS PECHCV AND PECHFV 
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defined, one for data that exists on only one vertical 

level, and one for data that exists on numerous vertical, 

or K, levels. These subscript functions use the row 

pointers, pointers to the variables in a row which are 

known as offsets, and the vertical level, k, for the second 

type of subscript function. 

The two subscript functions, which take the form 

of FORTRAN statement functions, are defined below. 

ISUBS (J,IOFF) = J*NVAR+IOFF IIII.7] 

for data on only one vertical level and 

ISUBK(J,IOFF,K) = J*WVAR+I0PF+(K-1)*IMAX [III. 8] 

for data on many vertical levels. In these definitions, 

J = Pointer to the row being referenced in the 
buffer array 

lOFF = Offset of the particular variable desired 
K = Vertical level desired 
WVAR = Total number of words per row record 

IMAX = Total length of a row {range on I) 

It is assumed that the pointers to the AO array {JMN, JN, 

JPN and JIN) range between 0 and 3 and that the pointers 

to the AN array (JMNP and JNP) range between 0 and 1. These 
subscript functions then point to the 1=0 element for 
ISUBS and the 1=0 element on level K for ISUBK. 

A simple example of the use of the subscript func- 
tions is called for at this point. Assume that it is desired 
3u 

to compute ^ on line J, level K; then the FORTRAN code for 
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doing this would be as follows: 

DO 10 K=1,KMAX 
IU=ISUBK(JN,DN,K) 

DO 10 I=2,IMAXM 

DU(I,K)=(AO{IU+I+l) -AO(IU+I-l) )/TDX 
10 CONTINUE 

The mechanics for accessing data in the storage 
buffers AO and AN have now been constructed. The data access 
mechanism together with the buffer structure shown in Figure 
III-2 allows a buffered computation to be performed; that 
is, simultaneous computations and I/O to riie files lOLD 
and INEW. The next, and final, step is to put these ideas 
together to generaue the data flow for one complete time 
step. The data flow shown assumes that free slip, insulated 
wall boundary conditions are used. Figure III- 3 shows the 
data flow for one sweep of the data file for the forecast 
models PECHCV and PECHFV. 

Finally, Tables III-l and III-2 define the offsets 
which are used by the two forecast models PECHCV and PECHFV 
in the pointer scheme to access data in the buffer arrays AO 
and AN. 
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Define Pointers JMN, JN, JPN, JIN 
for AO and JMNP, JNP for AN 


Bring in First 3 Rows from lOLD 
to AO, Lines JMN, JN, JPN 


Compute Lower Boundary Conditions 


Copy Row JMN in AO to Row JMNP 
in AN 


Start Buffer of Row JMNP in AN 
to INEW 


Compute Left-Right Boundary 
Condititons 


Compute Upper Boudna3ry Conditions 
If Last Time Through 


Copy Row JN in AO to Row JNP in 
AN, Update Below 


Compute Diagnostics for Row JN 
in AO 


Compute on Row JN in AO Storing 
Updated Data in Row JNP in AN 


Check Buffer of Row JIN in AO from 
lOLD, Except on Last Time Through 


Check Buffer of Row JMNP in AN to 
INEW 


FIGURE III-3: 


DATA FLOW FOR ONE SWEEP OF 
THE DATA FILE - MODELS PECHCV 
AMD PECHFV 
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Update Pointers JMN, JN, JPN, JIN 
for AO and JMNP, JNP for AN 


Test for End 

of Grid Sweep 


^ yes 


Write Out What Was Row JNP in AN 
to INEW 


'Write Out What Was Row JPN in AO 
to INEW (Boundary Row) 


Rewind lOLD and INEW 


Flip Names of lOLD and INEW 


FIGURE III-3: (Continued) 
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TABLE III-l: OFFSETS OF VARIABLES IN ARRAYS 

AO AND AN FOR PECHCV 


DESCRIPTION 

NAME 

LEVELS 

NO. 

Initial terrain pressure, 

PTI 

1 

0 

ir reduction to sea level 

PIR 

1 

63 

Terrain geopotential 

GZ 

1 

126 

Albedo and land-sea-ice 
indicator (packed) 

ALSI 

1 

189 

Sea surface temperature 

SST 

1 

252 

Tendency truncator 

RC 

1 

315 

Map factor, m 

MF 

1 

378 

Sine of latitude, sintj) 

SL 

1 

441 

Geopotential , * 

PHI 

5 

504 

Vertical velocity, w 

WV 

4 

819 

Precipitation 

PRE 

1 

1071 

Sea surface vapor pressure 

QSS 

1 

1134 

Saturation vapor pressure 

QS 

3 

1197 

Heating rates 

HT 

5 

1386 

X direction pressure gradient 
force, n and n-1 time levels 

PXN 

5 

1701 

Y direction pressure gradient 
force, n and n-1 time levels 

PYN 

5 

2016 

U wind, n and n-1 time 
levels 

DN 

5 

2331 

V wind, n and n-1 time 
levels 

VN 

5 

2646 


TABLE II I- 1: 


(Continued) 


DESCRIPTION 


NAME 

LEVELS 

i 

o 

Temperature, n and n 
levels 

“1 time 

TN 

5 

2961 

Moisture, n and n-1 
levels 

time 

QN 

3 

3276 

Terrain pressure, it, 
n time levels 

n+1 and 

PTNP 

1 

3465 

Terrain pressure, ir, 
n-1 time levels 

n and 

PTN 

1 

3528 



NVAR = 


3591 
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T2^LE III-2; OFFSETS OF VARIABLES IN ARRAYS 
AO AND AN FOR PECHFV 


DESCRIPTION 

NAME 

LEVELS 

NO. 

Initial terrain pressure, to, 
and t reduction to sea level 
(packed) 

PTIR 

1 

0 

Map factor, m 

MF 

1 

63 

Geopotential , ® 

PHI 

10 

126 

U wind, n and n-1 time levels 

UN 

10 

756 

V wind, n and n-1 time levels 

VN 

10 

1386 

Temperature, n and n-1 time 
levels 

TN 

10 

2016 

Moisture, n and n-1 time 
levels 

QN 

6 

2646 

Terrain pressure, ir, n+1 and 
n time levels 

PTNP 

1 

3024 

Terrain pressure, it, n and 
n-1 time levels 

PTN 

1 

3087 

Terrain geopotential 

GZ 

1 

3150 

Albedo and land-sea-ice 
indicator (packed) 

ALSI 

1 

3213 

Sea surface temperature 

SST 

1 

3276 

Sea surface vapor pressure 

QSS 

1 

3339 

Tendency truncator 

RC 

1 

3402 

Sine of latitude, sintf) 

SL 

1 

3465 

Precipitation 

PRE 

1 

3528 
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b. Data Storage - PEFHCV and PEFHFV 

The magnitude of the data file that must be proc- 
essed for the fine horizontal resolution forecast models 
is considerably larger than for the coarse horizontal res- 
olution models. This rather severely restricts the type of 
data flow scheme that can be constructed, since the problem 
becomes one of fitting the computation into the machine (a 
CYBER 76 was used as the target machine) . 

In order to fit the computation in the machine, a 
split row was devised, A split row contains all of the 
variables for a vertical slice, but they are grouped into 
two segments. The lower portion of the split row contains 
the variables that are required on the three rows J-1, J and 
J+1 during computation on row J. The upper portion of the 
split row contains variables that are only needed on row J. 
This allows smaller central memory buffers since there need 
be only one full row and two partial rows in central memory 
simultaneously. 

A diagram of the necessary central memory and ECS or 
LCM storage buffers for this scheme to work for the model 
PEFHCV is given in Figure III-4. There are three rows in 
central memory, a partial j-1 row, AJM, a full J row, AJ, 
and a partial J+1 row, AJP. There are also two buffer arrays 
IBUF and OBUF in central memory used for buffering data from 
and to files lOLD and INEW, respectively. For ECS or LCM 
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Central Memory 


AJM/ AJ and AJP arrays 



FIGURE III-4: STORAGE BUB’PERS REQUIRED 

FOR ROW STORAGE SCHEME 
IN MODEL PEFHCV 



III-22 







there is a buffer array, AIN, for three full rows. This 
buffer array is a circular buffer and has pointers JLM, 

JL and JLP pointing to rows J-1, J and J+1, respectively. 

Similarly a diagram of the necessary central 
memory and ECS or LCM storage buffers for this scheme to 
work for the model PEFHFV is given in Figure III-5. There 
are three rows in central memory, a partial J-1 row, AJM, a 
full J row, AJ, and a partial J+1 row, AJP. IN ECS or LCM 
there is a buffer array, AIN, for three full rovvs. This 
buffer array is a circular buffer and has pointers; JLM, JL 
and JLP pointing to rows J-1, J and J+1, respectively. Due 
to the way that the CDC system software works (buffer 
transfers into LCM are not allowed) , access to the files 
lOLD and INEW is made through the central memory array A-I 
with a block transfer either to or from ECS or LCM also 
being required. 

One more construct is needed before a given piece 
of data can be accessed in the arrays AJM, AJ and AJP. 

This is a method of pointing to the particular variable 
that is required. The offset of the variable can be used 
to point to single vertical level data since there are no 
row pointers to the central memory arrays. However, for 
computational convenience, a subscript function has been 
defined for data that exists on more than one vertical level, 
or K level. This STibscript function, which takes the form 




Central Memory 


from lOLD — ^ 
to INEW-^ — 


partial JP row 


full J row 


partial JM row 


AJM, AJ and AJP arrays 


Extended Memory or LCM 



AIN array 

Pointers JLM, JL and JLP 


FIGURE III-5; STORAGE BUFFERS REQUIRED 
FOR ROW STOBiAGE SCHEME 
IN MODEL PEFHFV 
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[III. 9] 


■s 


of a FORTRAN statement function, is defined below: 


ISUBK{IOFF,K) = IOFF+ (K-1) *IMAX 

where 

lOFF = Offset of the particular variable desired 

K = Vertical level desired 

IMAX = Total length of a row (range on I) . 

For single level data the offset of the variable points to 
the 1=0 element, and for multiple vertical level data the 
subscript function ISUBK points to the 1=0 element for 
level K. An example of the use of the subscript function 
was given in the previous section. 

The mechanics for accessing data in the central 
memory arrays AJM, AJ and AJP have now been constructed. 

The data access mechanism together with the buffer structures 
shown in Figures II1-4 and III-5 allow the computation to 
sweep the data base on file lOLD while constructing an 
updated data base on file INEW. Figures III"6 and III-7 
combine these ideas and present the data flow for one com- 
plete sweep of the data file in the forecast models PEPHCV 
and PEFHFV, respectively. The data flow shown assumes that 
free slip, insulated wall boundary conditions are used. 

Finally, Table III-3 defines the offsets which are 
used by the two forecast models to access data in the buffer 
arrays AJM, AJ and AJP, 
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FIGURE III-6: DATA PLOW FOR ONE SfJEEP OP THE 

DATA FILE - MODEL PEFHCV 
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FIGURE III-6: (Continued) 
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FIGURE III-6: 


(Continued) 






Define pointers JLM, JL and JLP 
to array AIN in LCM 


Buffer in first 3 rows to AJ from 
lOLD, block transfer to rows JLM, 
JL and JLP of AIN in LCM 


Block transfer first 3 rows from 
AIN in LCM to AJM, AJ and AJP 
in SCM 


Compute lower boundary conditions 


Block transfer AJM in SCM to row 
JLM of AIN in LCM 


Block transfer row JLM of AIN in 
LCM to AJ in SCM 


Buffer out first row from AJ to 
INEW 


Initialize line counter LINE=2 


Block transfer rows JLM, JL and 
JLP of AIN in LCM t.T AJM, AJ and 
AJP in SCM 


Compute E-W boundary conditions 


LINE . NE . JMAXM 


Compute upper boundary conditions 


Compute on row J in AJ 


FIGURE III-7: DATA FLOW FOR ONE SWEEP OF 

THE DATA FILE - MODEL PEFHFV 
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FIGURE III-7: (Continued) 
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TABLE II 1-3 : 


OFFSETS OF VARIABLES IN ARRAYS 

AJM^ AJ AND AJP FOR PEFHCV AND PEFHFV 


DESCRIPTION 

NAME 

5 LAYER 
Level No. 

10 LAYER 
Level No. 

Initial terrain pressure, 
TTof and TT reduction to 
sea level (packed) 

PTIR 

1 

0 

1 

0 

Map factor , m 

MF 

1 

187 

1 

187 

Geopotential heights , $ 

PHI 

5 

374 

10 

374 

U wind, n and n-1 time levels 

UN 

5 

1309 

10 

2244 

V wind, n and n-1 time levels 

VN 

5 

2244 

10 

4114 

Temperature, n and n-1 
time levels 

TN 

5 

3179 

10 

5984 

Moisture, n and n-1 
time levels 

QN 

3 

4114 

6 

7854 

Terrain pressure, tr, n+1 
and n time levels 

PTNP 

1 

4675 

1 

8976 

Terrain pressure, v , 
n and n-1 time levels 

PTN 

1 

4862 

1 

9163 

Terrain geopotential 

GZ 

1 

5049 

1 

9350 

Albedo and land-sea-ice 
indicator (packed) 

ALSI 

1 

5236 

1 

9537 

Sea surface temperature 
and sea surface vapor 
pressure (packed) 

SSTQSS 

1 

5423 

1 

9724 

Tendency truncator 

RC 

1 

5610 

1 

9911 

Sine of latitude, sin4> 

SL 

1 

5797 

1 

10,098 

Precipitation 

PRE 

1 

5984 

1 

10,285 
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TABLE III- 3 ; OFFSETS OF VARIABLES (Continued) 


DESCRIPTION 

NAME 

5 LAYER 
Level No. 

10 LAYER 
Layer No. 

Heating Rates 

HT 

3 

6171 

3 

10,472 

Vertical velocities^ w 

W 

4 

6732 

9 

11,033 

X direction pressure 
gradient force ^ n and 
n-1 time levels 

PXN 

5 

7480 

10 

12,716 

Y direction pressure 
gradient force, n and 
n-1 time levels 

PYN 

5 

8415 

10 

14,586 


NVARS = 


5049 


9,350 


NVARL = 


9350 


16,456 


NOTE: 


IMAX = 187 


B. Code Structure - PECHCV and PECHFV 


This section describes the program structure of the 
forecast models PECHCV and PECHFV. It also defines the 
various files used by the programs, defines the common 
blocks that are used and describes each of the subprograms 
that constitute the models . 

The codes are written in CDC FORTRAN, which is a dia- 
lect of FORTRAN IV, for the CDC FORTRAN Extended Compiler 
(0PT=2 level) . The two programs, PECHCV and PECHFV, have 
been designed to use rotating mass storage communicating 
with central memory through two I/O channels. 

The forecast models are fairly large codes. Since the 
various sections of the codes — initialization, integration 
and output — are each fairly large, the codes have been 
overlaid to more efficiently use central memory. The sec- 
tions that follow describe the main divisions of the fore- 
cast models by overlays. Section 1 describes the control 
overlay, gives a block diagram of the overlay structure, 
defines the I/O file structure and defines the common blocks 
that are global to all overlays. 

Section 2 describes the initialization overlay, defines 
the common blocks local to this overlay, and describes each 
of the subprograms that are used by this overlay. Similarly, 
Sections 3 and 4 describe the integration and output over- 
lays, respectively. 
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1. Main Overlay - Programs PECHCV and PECHFV. 

The main overlay of the forecast model programs is 
very simple, since it only performs the function of declar- 
ing three global common blocks, defining various control 
variables and '■ailing the actual working overlays. The 
overlay structure of the programs is shown in Figure III-8. 
Only one level of overlay is used for programs PECHCV and 
PECHFV. 

The various files used by the programs are defined in 
Table III-4. Tables III-5 and III-6 define the variables 
that appear in the global common blocks CNTRL, and INDEX, 
respectively. Common block BUFPOINT consists of the offsets 
for the variables and contains the integer variables listed 
in Tables III-l and III-2 for PECHCV and PECHFV, respectively. 

a. Subroutine BUFERR 

Subroutine BUFERR prints a diagnostic message 
when an error occurs on a BUFFER (IN or OUT) operation that 
tells where the error occurred. 
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TABLE III-4: FILE STRUCTURE OF FORECAST 

MODELS PECHCV AND PECHPV 


OUTPUT Used to write various diagnostic 

messages 

TAPEl Input data file for initialization 

overlay; generated by analysis 
programs 

TAPE2 Random access file used by 

initialization and output 
overlays 

TAPE3 Restart data tape; contains a \ 

header record and contents of i 

current lOLD file ; 

TAPE4 Output data file; contains fore- | 

cast fields to be selectively \ 

plotted by graphics program and for 
use by the analysis programs 

TAPES {lOLD, INEW) integration data 

TAPES (lOLD, INEW) integration data | 
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TABLE 

lOLD 

INEW 

ISW 

ID (2) 

LIMI 

LIM 

ITAU 

ITSPHR 

NN 

NNOUT 

LCO 

L12 

DAY 

HR 

DT 

TDT 

MATS 


III-5: DEFINITION OF VARIABLES IN COMMON BLOCK 

/CNTRL/ FOR PECHCV AND PECHFV 


File name of integration data 
(TAPES or TAPES) 

File name of integration data 
(TAPES or TAPES) 

=1/ SSWl is "on", signifying normal 
start with data from TAPEl 
=2, SSWl is "off", signifying restart 
with data from TAPES 

1st word, forecast identifier (date-time group) 
2nd word, forecast data field identifier 

Starting cycle 

Number of integration and output cycles 
(2 X number of days + 1) 

Forecast time indicator (hours) 

= 0 during initialization and first output cycle 
= 12*LC0 for normal integration cycle 

Number of integration time steps per hour 

Time step counter, reset after each output 
cycle, also NN=0 signifies a Matsuno time step 

Number of time steps between normal 
output cycles, usually equal to 
12 X ITSPHR 

Output cycle counter, updated by 1 each 
output cycle 

Indicator set at 12~hour output cycle, 
used to zero rainfall 

Number of days from January 1 
of the forecast 

Hour of the day, GMT 

Time step. At 

Twice the time step, 2At 

Logical variable used during Matsuno time step 
= .FALSE, during first half 
= .TRUE, during second half 
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TABLE III-6: DEFINITION OP VARIABLES IN COMMON BLOCK 

/INDEX/ FOR PECHCV AND PECHPV 


ICL 

ICU 


IMAXM 

IMAX 

JMAXM 

JMAX 

KMAXQ 

KMAXM 

KMAX 

LINE 

NFLD 

NVAR 


Start of computation on data row (=2) 

End of computation on data row 
(=IMAX-1) 

Length of data row -1 

Length of data row 

Number of data rows -1 

Number of data rows 

Number of moisture levels 

Number of vertical levels -1 

Number of vertical levels 

Number of row being worked on 

Size of a data record on TAPE2 (63 x 63) 

Size of a data record on TAPES and TAPE 9 


2 . 


Initialization Overlay - Program INITPE 


The function of the initialization overlay is to 
take as input data the fields generated by the analysis 
programs and to generate an initial state for the inte- 
gration overlay which actually performs the forecast. That 
is, the interpolation and synthesis of data fields to the 
PEM sigma surfaces using data produced by the analysis 
programs as described in Section II-5 are performed by 
the initialization overlay. 

Program INITPE extracts data generated by the analysis 
programs from the sequential file TAPEl. Using these data, 
it interpolates and synthesizes the data fields on the sigma 
surfaces and writes these new fields to the random file 
TAPE2. They are then sorted into vertical slice format on 
file lOLD for use by the integration overlay. 

Figure III-9 shows the structure of the initialization 
overlay, giving the names of all the subprograms which 
constitute this overlay. The following sections describe 
these routines. 

The labeled common block, INCON, is local to this 
overlay and its variables are described in Table III-7. 

Table III-8 defines the data fields on TAPEl which are 
produced by the analysis programs. Table III-9 defines the 
fields on TAPE2 which are produced both by the initiali- 
zation process to define the initial state of the forecast 
and by the output overlay to present the results of the 
forecast. 
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TABLE III- 

PLVL(13) 

PLVLN (13) 
SIG(5)or(10) 

ZLN(5)or(10) 
SLN (5)or (10) 
G 

GRDC 

CTOK 

AE 

BE 


: DEFINITION OF VARIABLES IN COMMON BLOCK 

/INCON/ FOR PECHCV AND PECHFV 


Pressure levels of input analysis data: 

1000, 950, 900, 850, 700, 500, 400, 300, 

250, 200, 150, 100, 50 

£n [PLVL(K)1, K = 1,KMAX 

Sigma levels 0.9, 0.7, 0.5, 0.3, 0.1 
for PECHCV and .95, .85, .75, .65, .55, 

.45, .35, .25, .15, .05 for PECHFV 

2-n [SIG(K)], K = 1,KMAX 

R*iln(0.9) for K=l, R*2n [SIG (K) /SIG (K-1) ] ; K=2 

2 

Gravitational acceleration = 9.80616 m/sec 

Constant used in calculation of sin (latitude) 
(31.205)^ 

Conversion constant relating degrees 
Celsius to degrees Kelvin =273.16 

Constant used in calculation of saturation 
vapor pressure =21.656 

Constant used in calculation of saturation 
vapor pressure = 5418 


,KMAX 


TABLE III-8 


CONTENTS OP THE INPUT 

DATA FILE - TAP El 

CREATED BY THE ANALYSIS PROGRAMS 


Record 


Data Field 


1 

Geopotential heights 

of the 

1,000 mb 

surface 

2 

It 

11 


950 

II 

3 

IT 

It 


900 

II 

4 

It 

It 


850 

tt 

5 

tl 

II 


700 

II 

6 

II 

II 


500 

11 

7 

II 

It 


400 

IT 

8 

It 

II 


300 

II 

9 

11 

It 


250 

Tl 

10 

tl 

11 


200 

IT 

11 

IT 

II 


150 

II 

12 

It 

II 


100 

11 

13 

Tt 

ir 


50 

IT 

14 

Temperatures 

on the 1 

,000 mb 

surface 


15 

II 

II 

950 

If 


16 

tt 

tt 

900 

II 


17 

JT 

11 

850 

tt 


18 

It 

If 

700 

U 


19 

11 

II 

500 

II 


20 

II 

IT 

400 

11 


21 

It 

11 

300 

It 


22 

II 

II 

250 

II 


23 

11 

II 

200 

II 


24 

n 

IT 

150 

11 


25 

u 

II 

100 

n 


26 

II 

11 

50 

IT 



27 

28 


Land-sea-ice indicator 
Surface pressure, Pg 
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TABLE III -St INPUT DATA FILE - TAPEl (Continued) 


Record 


Data Field 


29 Terrain height 

30 Sea surface temperature 

31 Albedo 


32 

U 

wind 

on the 1,000 

mb surface 

33 

V 

11 

1,000 

II 

34 

U 

11 

950 

M 

35 

V 

It 

950 

II 

36 

U 

11 

900 

fl 

37 

V 

11 

900 

11 

38 

U 

11 

850 

tl 

39 

V 

K 

850 

Tl 

40 

U 

11 

700 

II 

41 

V 

II 

700 

II 

42 

U 

tl 

500 

11 

43 

V 

IT 

500 

11 

44 

U 

It 

400 

11 

45 

V 

ft 

400 

IT 

46 

U 

11 

300 

11 

47 

V 

II 

300 

tl 

48 

U 

tl 

250 

It 

49 

V 

II 

250 

II 

50 

u 

11 

200 

ri 

51 

V 

It 

200 

It 

52 

u 

It 

150 

It 

53 

V 

II 

150 

It 

54 

U 

It 

100 

11 

55 

V 

It 

100 

II 


NOTE; Geopotential heights are in meters; 

temperatures are in degrees centigrade; 
surface pressure is in mb; 
terrain height is in meters; 
winds are in meters per second. 

All fields are 63 x 63 arrays with a 
20 word ID prefixed. 
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TABLE 


Record 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 


III-9: CONTENTS OP THE RANDOM ACCESS FILE - 

TAPE2 FOR PECHCV AND PECHPV 10, SIGMA 
LEVEL VERSION OP THE FORECAST MODEL 


Data Field 


Analyzed geopotential heights 


of 1000 rab surface {meters) 
950 
900 
850 
700 
500 
400 
300 
250 
200 
150 
100 
50 


Analyzed temperatures of 1000 mb surface (=C) 
" 950 

” 900 

" 850 

•' 700 

" 500 

" 400 

" 300 

'• 250 

" 2 00 
" 150 

" 100 

" 50 


Analyzed surface pressure, Ps (mb) 
Terrain height, Z,p (meters) 

Sea surface temperature (°C) 
Land-sea-ice indicator 
Sin of latitude, sin<{> 

Map factor, m 
Tendency truncator 
Albedo 

Terrain geopotential, gZrp 


TABLE III-9: RANDOM ACCESS FILE ^ TAPE2 

(Continued) 


Record Data Field 


41 

Temperature on a = 0.95 surface 

(°K) 

42 

ir 

0.85 " 



43 

II 

0.75 •' 



44 

It 

0.65 '• 



45 

11 

0.55 



46 

tf 

0.45 



47 

11 

0 - 35 " 



48 

11 

0,25 



49 

(1 

0.15 " 



50 

11 

0.05 " 



51 

Vapor pressures on a = 0,95 

surface (mb) 

52 

II 

0.85 

It 


53 

ft 

0.75 

If 


54 

11 

0.65 

II 


55 

If 

0.55 

tf 


56 

II 

0.45 

II 


57 

Terrain 

pressure, r (mb) 



58 

Pressure 

difference, Ps—ir 

(mb) 


59 

An (tt) 




60 

Forecast 

Ps (mb) 



61 

Geopotentials of o =0,95 surface 

(g« meters 

62 

n 

0.85 

II 


63 

ir 

0.75 

tt 


64 

It 

0.65 

II 


65 

ir 

0.55 

n 


66 

11 

0.45 

II 


67 

It 

0.35 

u 


68 

t( 

0.25 

II 


69 

tt 

0.15 

Tt 


70 


0,05 

II 


71 

U winds 

on 950 mb surface 


(m/sec) 

72 

tt 

850 



73 


700 



74 

tl 

500 



75 


400 



76 

Iff 

300 



77 

m 

250 



78 

II 

200 



79 


150 



80 


100 




rr:- 
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I 



TABLE III- 9: 


RANDOM ACCESS FILE - TAPE2 
(Continued) 


Record 


Data Field 


81 

V 

winds 

on 950 

mb surface 

(ra/sec) 


82 


It 

850 





83 


II 

700 





84 


tl 

500 





85 


tl 

400 





86 


It 

300 





87 


II 

250 





88 


tl 

200 





89 


It 

150 





90 


II 

100 





91 

U 

\7inds 

on a = 

0.95 

surface 

(m/sec) 

92 


II 


0,85 

It 



93 


tl 


0.75 

Tl 



94 


II 


0.65 

II 



95 


tl 


0,55 

II 



96 


11 


0.45 

11 



97 


tt 


0.35 

II 



98 


II 


0.25 

11 



99 


11 


0.15 

It 



100 


tl 


0.05 

II 



101 

V 

winds 

on 0 = 

0.95 

surface 

(m/sec) 


102 


Tl 


0.85 

ri 



103 


II 


0,75 

tl 



104 


II 


0.65 

ir 



105 


ir 


0.55 

Cl 



106 


II 


0.45 

IT 



107 


II 


0.35 

It 



108 


11 


0.25 

It 



109 


tt 


0.15 

II 



no 


II 


0.05 

It 



111 Forecast 

geopotential 

height of 1, 

000 mb surf 

112 

H 


II 



950 

11 

113 

II 


II 



900 

II 

114 

It 


tl 



850 

II 

115 

II 


tl 



700 

tt 

116 

It 


II 



500 

tl 

117 

tl 


It 



400 

It 

118 

II 


II 



300 

A 

119 

II 


U 



250 

II 

120 

n 


11 



200 

It 

121 

tl 


It 



150 

tl 

122 

II 


IT 



100 

II 


(meters) 
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TABLE III-9: RANDOM ACCESS FILE - TAPE2 

(Continued) 


Record 

123 

124 

125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 

161 

162 

163 

164 

165 


Data Field 


Forecast temperatures on 

Tl 

It 

ft 

V 

It 

II 

u 

It 

It 

It 

It 


1,000 mb surface (®C) 
950 '• 

900 

850 " 

700 

500 

400 

300 

250 

200 " 

150 *' 

100 


Forecast U winds on 1,000 


It 

11 

It 




950 

900 

850 

700 

500 

400 

300 

250 

200 

150 

100 


mb 


surface 

II 


It 

II 


II 


IT 

n 

II 


(m/sec) 


Forecast V winds 

on 1,000 mb 

surface 

Tl 

950 

It 

11 

900 

11 

11 

850 

11 

It 

700 

It 


500 

tt 

II 

400 

II 

II 

300 

11 

II 

250 

11 

IT 

200 

tl 

n 

150 

tt 

It 

100 

II 


(m/sec) 


Sea surface saturation vapor pressure 
Forecast precipitation (cm) 

Vapor pressure on 1,000 mb surface 
" 900 “ 

II yOO II 

500 " 

'• 300 ” 


(mb ) , 
(mb) 
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a. 


Subroutine SETPHIS 


Subroutine SETPHIS computes the initial geopoten- 
tials of the sigma surfaces. Input to the routine consists 
of the temperatures on the sigma surfaces and the moistures 
on the lower sigma surfaces. The geopotentials are computed 
using the form of the hydrostatic equation given by Equa- 
tions [11.107] and [11.108] given in Section II-D-1. 

Finally, the computed geopotentials are stored on the 
random access file TAPE2 for use in the forecast. 

b. Subroutine MOIST 

Subroutine MOIST computes the relative humidity 
on a pressure surface using the procedure based on the rel- 
ative vorticity as described in Section II-E-6. Input to 
the routine consists of the geopotential heights of the 
pressure surface. The routine then computes the relative 
vorticity on the pressure surface and, using the empirical 
relation given by Equation [11.133], computes the relative 
humidity on the pressure surface. 

Program INITPE calls this subroutine five times 
to generate moistures on the 1,000, 900, 700, 500 and 
300 mb surfaces. Data on these pressure surfaces is then 
used to obtain initial moisture data on the lower sigma 
surfaces through logarithmic pressure interpolation. 


These data fields are stored on the random access file 
TAPE2 for use in the forecast. 

c. Subroutine WIMP 

Subroutine WIND computes the geostrophic wind 
components at 50 mb for use in the initialization of the 
winds on the model sigma surfaces. Input to the routine 
consists of the geopotential heights of the 50 mb pressure 
surface. 


d. Subroutine TREAD 

Subroutine TREAD acts as a generalized input data 
routine to obtain analysis data from TAPEl. A storage array 
address and the data field ID are input to TREAD from the 
calling routine. TREAD then searches a table (array ITP) 
to find the location of the desired data field and reads 
the data field from TAPEl into the storage array. Various 
data field orderings on TAPEl may be accommodated simply 
by changing the information stored in array ITP. 

e. Subroutine DATAIQ 

Subroutine DATAIO is a generalized routine for 
writing and reading data fields to and from the random 
access file TAPE2. Two entry points exist for this routine: 
DATAOUT, which uses the system routine WRITMS for writing to 
the file, and DATAIN, which uses the system routine READMS 
for reading from the file. 
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f . Subroutine RSTART 


Subroutine RSTART reads the restart data tape, 
TAPES ^ when a forecast run is being restarted and initial- 
izes the forecast model in order to continue the forecast. 
Restarts are controlled by sense switch 1, If SSWl is 
"ON", then a normal run is being made; if SSWl is "OFF", 
then a restart run is being made. 

A restart data tape is written during each output 
cycle so if the computer goes down or if for some reason a 
forecast is to be continued at a later time, the forecast 
can be resumed from the point at which the last restart tape 
was written. Table III-IO gives the format of the restart 
data tape TAPE3. 

Thus subroutine RSTART reads in the header record 
from TAPES and initializes the required control variables. 
Next the routine reads in the forecast data records and 
recreates file lOLD. The forecast run will now proceed from 
this point in the forecast. 

g. Subroutine SORTIT 

Since the output of the analysis programs and the 
initialization routines are data on horizontal planes and 
the integration overlay performs computations on data 
arranged in vertical slices , an interface routine to pro- 
perly structure the data is required. Subroutine SORTIT 


TABLE III-IO: RESTART DATA TAPE (TAPES) FORMAT 


1) Header record (twenty words) - 

contains information to identify the forecast 


1 

ID(1) 

2 

ITAU 

3 

lOUT 

4 

LCO 

5 

L12 

6-10 

Unused 

11 

DAY 

12 

HR 

13-20 

Unused 


2) Forecast data records - 

these contain the contents of file lOLD at ITAU 
(JMAX records in vertical slice format) 


takes the initial data in horizontal plane form which is 
stored on the random access file TAPE2 and constructs the 
forecast data file lOLD in vertical slice form for use by 
the integration overlay. 

Due to core storage limitations, more than one 
pass must be made through the file lOLD to completely 
initialize the forecast. Each pass consists of reading a 
maximirm of ten planes of data from TAPE2 and then sweeping 
file lOLD from row one to row JMAX, writing out the initial 


data. 


3. 


Integration Overlay - Program INTGPE 


The function of the integration overlay is to numer- 
ically integrate the primitive equations; that is^ to 
step the data fields forward in time to produce a forecast. 
Input to this overlay consists initially of the data field 
specification produced by the initialization overlay and/ 
subsequently/ the current state of the forecast contained 
on file lOLD. After initialization, a forecast consists of 
cycling through the integration and output overlays until 
the desired forecast time has been reached. 

Figure III-IO shows the structure of the integration 
overlay, giving the names of all the subprograms which 
constitute this overlay. The function of each of these 
routines is described in the sections that follow. 

Table III-ll defines the variables that appear in the 
common blocks that are local to the integration overlay for 
programs PECHCV and PECHPV. 

a. Subroutine DEFINE 

Subroutine DEFINE has the simple function of 
defining various constants that are used in the integration 
overlay and zeroing the buffer arrays of this overlay. The 
subroutine is called once per integration-output cycle upon 
entering the integration overlay. 
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TABLE III-ll: DEFINITION OF VARIABLES IN 

COMMON BLOCKS LOCAL TO THE 
PROGRAMS PECHCV AND PECHFV 


Common 


Common 


Common 


Block /BUFF/ 

A0{NVAR,4) Circular input buffer for rows 

JMN, JN, JPN and JIN 

AN(NVAR/2) Circular output buffer for rows 

JMNP and JNP 


Block /TEND/ 

UT(63,KMAX) 

VT(63,KMAX) 

TT (63,KMAX) 

QT(63,KI1AXQ) 

PIT (63) 


, tendencies in iru for all 

levels on row JN 

3ttv 
3t 

levels on row JN 

3irT 
3t 

levels on row JN 


f tendencies in itv for all 


, tendencies in itt for all 


, tendencies in 'jrq for all 


3t 

levels on row JN 

, tendency in terrain pressure, it, 
on row JN 


Block /LINEAR/ 

DPIDX (63,KMAX) Kurihara modification pressure 

gradient terms in x direction for all 
levels on row JN 

DPIDY (63,KMAX) Kurihara modification pressure 

gradient terms in y direction for all 
levels on row JN 

W(63,KMAXM) Vertical velocities for all levels 

of row JN 



TABLE III-ll: Definition of Variables (Continued) 

Coimnon Block /LOCAL/ 

UPM(6 3) 

VPM(63, 3) 

PILN(63, 3) 

PIKAP (63) 

TCOR(63; 3,KMAX) temperatures corrected for vapor 

pressure at all levels on rows JMN, 

JN and JPN 

NOTE ; TCOR is equivalenced to the 
following variables: 

Arrm^ 

AB(63,KMAX) (6a + 5 3) at all levels on row JN 

QC(63fKMAXQ) vapor pressure, q, at all levels on 

row JN 

QCS (63 ,KMAXQ) saturation vapor pressure, q , at all 

levels on row JN 

TC(63,KMAX) temperature at all levels on row JN 

PC(63,KMAXQ) pressure at moisture levels on 

row JN 

precipitation on row JN 


— at a particular level on row JN 

at a particular level on rows 
JMN, JN and JPN 

In(Tr) on rows JMN, JN and JPN 

on row JN 
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RF (63) 


TABLE III-ll: DEFINITION OF VARIABLES (Continued) 


Conation Block /CON/ 


D 

Ax = Ay = 381000 for coarse-mesh models and 
127000 meters for fine-*mesh models 

DSQ 

Ax^ = Ay^ 

TD 

2Ax = 2Ay 

OPD 

1/4 Ax 

DSIG 

Act = 0.2 for five— layer models and 0.1 for 
ten— layer models 

TDSIG 

2Act 

SIG(KMAX) 

Sigma levels 0.9, 0.7, 0.5, 0.3, 0.1 for five- 
layer models, and 0.95, 0.85, 0.75, 0.65, 0.5! 
0.45, 0.35, 0,25, 0.15, 0.05 for ten— layer 
models 

CT (KMAXM) 

In [SIG(K)/SIG(K“D] ; K = 2, KMAX 

SLN (KMAX) 

R*ln(0.9) or In (0.95) for K=l, R*CT(K) for 
K=2,KMAX 

S9* 

In (0.9) 

S7* 

In (0.7) 

S5* 

In (0.5) 

S3* 

In (0.3) 

SI* 

}n (0.1) 

AE 

Constant used in calculation of saturation 
vapor pressure = 21.656 

BE 

Constant used in calculation of saturation 
vapor pressure = 5418.0 

G 

2 

Gravitational acceleration = 9.80616 m/sec 



TABLE Ill-il: Definition of Variables (Continued) 


R 

XKAP 

FC 

CD 

DK 


CK(KI-IAX) 


Gas constant = 287.04 
K = Vcp = 0.2858 

Coriolis parameter = 20 ,- 1.45842x10"^ 

Drag coefficient = 0.0015 

6 2 2 

Diffusion coefficient = 1x10 in'^/sec'^ for 

5 

coarse-mesh models and 1x10 for fine-mesh 
models 

[1000/SIG(K)] **XKAP,K=1,KMAX 


* These five constants are appropriate to the five— layer 
models. They are replaced with S95=ln (0 . 95) , etc. 
for the ten— layer models. 


Note: Symbolic dimensions are given where appropriate 

since the size of the arrays varies with the forecast models 
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b. 


Subroutine LSTEPP 


Subroutine LSTEPP perforins the first half of a 
pressure-gradient-force averaged time stepr integrating 
forward one time step the pressure tendency, thermodynamic 
energy and moisture equations. The routine controls the 
data flow for one sweep of the data base and calls the 
routines necessary to compute diagnostics, boundary con- 
ditions, forcing functions and diabatic effects. Figure 
III-IO shows the controlling function of subroutine LSTEPP. 

P\^ure III-3, which appears in Section III-2-a 
shows the data flow for one sweep of the data base for models 
PECHCV and PECHPV. This figure concentrates on the data 
flow aspect of the routine and the computational portion is 
denoted by the box labeled "computation on row J". Figure 
III-ll is an expansion of this box for LSTEPP and shows the 
step-by-step computations and the subroutines that are 
called. 

c. Subroutine PITEND 

Subroutine PITEND computes the Eurihara pressure 
gradient terms for use in the thermodynamic energy equation, 
the pressure tendencies and the vertical velocities. The 
terrain pressure tendency is computed using Equation [11.104] 
and the vertical velocity is computed using Equation [II. 106] 
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Compute boundary condxtions (BNDRY) 


Compute diagnostics on row J {DIAG 
and DIAGL) 


Empty precipitation storage every 
12 hours 


Compute Air/ Ax, Ait/ Ay, pressure 
tendencies and vertical velocities 
(PITEND) 


Integrate terrain pressure 


Non-linear smoothing of terrain 
pressure (PISM) 


Compute tendencies for T and q 
(TQTEND) 


Integrate and timesmooth temperature 


Store q, q^, p and T into working 
arrays for diabatic computations 


Large scale condensation (LARCON) 


Compute radiation, evaporation and 
sensible heating on the hour (HEATING) 


Moist and dry convective adjustment 
(CONVADJ) 


Store modified q, q_ and T 


Compute Geopotentials (SETPHI) 


FIGURE III-ll: COMPUTATION FOR ONE TIME STEP IN LSTEPP 
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d. 


Siibroutine PISM 


Subroutine PISM performs the non-linear terrain 
pressure smoothing that is described in Section II-D-2. 

The non-linear smoothing operator given by Equation [11.118] 
is applied along an entire row immediately after the terrain 
pressure has been advanced in time. 

e. Subroutine TQTEND 

Subroutine TQTEND computes the tendencies for the 
temperature and vapor pressure using Equations [11.105] and 
[11.1091, respectively. 

The computation of the temperature tendency 
requires evaluation of the horizontal Arakiwa advective 
operator, the vertical flux and the adiabatic heating/ 
cooling term. Computation of the vapor pressure tendency 
involves the horizontal Arakawa advective operator and the 
vertical flux. 


f . Subroutine LARCON 

Subroutine LARCON computes the amount of precipi- 
tation resulting from the removal of excess moisture from 
super-saturated regions of the atmosphere and the temperature 
change due to the release of latent heat. The precipitation 
algorithm used ir. this computation is described in Section 
II-C-2. 
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g- 


Subroutine HEATING 


Subroutine HEATING perforins the diabetic heating 
computations which include shortwave and longwave radiation 
and sensible heating from the planetary boundary layer as 
described in Section II-C-1. Evaporation of moisture from 
the surface into the lowest sigma level is also computed. 

Since the radiation and boundary layer computations 
are complicated and time-consuming, they are only performed 
once every model-hour (normally every four time steps, when 
a fifteen-minute time step is used) . Temperature tendencies 
are then computed and saved for use over the next hour of 
integration. 


h. Function ALON 

Function ALON computes the longitude of a Northern 
Hemispheric polar stereographic (I,J) grid point in radians 
ranging from 0 - 2Tr clockwise from the Greenwich Meridian. 

i. Subroutine CONVADJ 

Subroutine CONVADJ performs the moist and dry 
convective adjustment as described in Section II-C-3. The 
moist convective adjustment consists of Arakawa's param- 
eterized cumulus convection which may result in precipita- 
tion and adjustment of the temperature and moisture in the 
lower three sigma layers for the five— layer models and in 


& 
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the lower six sigma layers for the— ten layer models. The 
dry convective adjustment process consists of checking the 
model atmosphere for static stability and modifying the 
temperature structure when necessary to achieve stability. 

j . Subroutine SBTPHI 

Subroutine SETPHI computes the new geopotentials 
after all the prognostic variables have been updated on the 
first half of the time step using Equations III.107] and 
[11.108] . 


k. Subroutine LSTEPV 

Subroutine LSTEPV performs the second half of a 
pressure-gradient— force averaged time step integrating the 
momentum equations forward one time step. The routine 
controls the data flow for one sweep of the data base and 
calls the routines necessary to compute diagnostics, boundary 
conditions and forcing functions. Figure III-IO shows the 
controlling function of subroutine LSTEPV. 

Figure III-3, which appears in Section III-2-a 
shows the data flow for one sweep of the data base for models 
PECHCV and PECHFV. This figure concentrates on the data 
flow aspect of the routine and the computational portion is 
denoted by the box labeled "computation on row J" . Figure 
III- 12 is an expansion of this box for LSTEPV and shows the 
step-by-step computation and the subroutines that are called. 
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Compute boundary conditions (BNDRY) 


Compute diagnostics on Row J 
(DIAGLP) 


Compute time— averaged pressure- 
gradient— force terras (PGPTRM) 


Compute tendencies for momentura 
equations (UVTEND) 


Modify lowest level tendencies for 
friction (DRAG) 


Integrate U wind 


Integrate V wind 


Move terrain pressure in central 
memory buffer 


Print diagnostics (DIAGP and DIAGPL) 


FIGURE III-12: COMPUTATION FOR ONE TIME 

STEP IN LSTEPV 
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1 . 


Subroutine PGPTRM 


Subroutine PGFTRM computes the time-averaged pressure 
gradient-force terras as given by Equation [II. Ill] for use 
in the momentum equations. The Kurihara modification is 
used to evaluate the terms on local pressure surfaces rather 
tiian directly on tlie sigma surfaces. 

The routine also stores the vertical velocities 
into a working array for use in subsequent computations. 

m. Subroutine UVTEMD 

Subroutine UVTEND computes the tendencies for the 
momentum equations using Equations [11.94] and [11.95]. 

The computation of the tendencies for the momentum 
equations involves the horizontal Arakawa advective operator ^ 
thf vertical flux, the time-averaged pressure-gi"adient— force 
term, the Coriolis force and the diffusion term. 

n. Subroutine DRAG 

Subroutine DRAG modifies the u and v momentum ten- 
dencies in the lowest layer to account for frictional dissi- 
pation in the gross boundary layer. The friction terms are 
evaluated using Equations [11.98] and [11.99], 
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o. 


Subroutine BNDRY 


Subroutine BNDRY computes the horizontal boundary 
conditions for the hemispheric polar stereographic grid. 

As explained in Section II-B, the horizontal boundaries con- 
sist of rigid impermeable vertical walls placed on the grid 
on the next to outermost row of grid points. The boundary 
conditions computed are given by Equation [II. 9]. 

Three entry points exist for this subroutine; 

LBNDRY, EWBNDRY and UBNDRY. LBNDRY computes values for the 
first row of the grid and is called once for each sweep of 
the data base. EWBNDRY computes values for the end points 
of each row of the grid and is required for each row as the 
grid is swept. UBNDRY computes values for the last row of 
the grid and is called once for each sweep of the data base. 

p. Subroutine DIAGN 

Subroutine DIAGN computes and prints various layer 
mean diagnostic quantities in order to monitor the forecast. 
Table III-12 defines the diagnostic quantities computed by 
DIAGN. In addition to the mean values computed by the routine, 
vertical column values are saved and printed for one arbi- 
trary point. 

There are three entry points in DIAGN: IDIAG, which 

initializes the diagnostic computations; DIAG, which computes 
diagnostics on a line-by-line basis; and DIAGP, which prints 
the layer mean diagnostics at the end of the time step. 
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TABLE III-12: DIAGNOSTIC QUANTITIES COMPUTED 

BY SUBROUTINE DIAGN 

□ _ 1 2 

^ “ N i "«ean of the u wind 

^ ifj -taycr moan of the v wind 

“ ifj mean of the temperature 

^ ifj of the goopotential 

Hr = — y ijm 

k N «^i,j,k/ layer mean of the diabetic heating 

^ ifj mean of tlie vortical velocity 

^ ifj mean of the vapor pressure 

"'^k - N Pi,j,k (“f,j,k + ^Lj,k>/<2'R*T.^j j^), 
layer mean of kinetic energy 

MvsQ^ = i j; 

i,j ‘“i+i,j,k - ■>• 

layer mean of the square of the divergence 

voiisQ. = i 2 r fv 

i,j i+l»j.k '^i-l,j,k ■ "i,j+i,k “i,j-i,i 


layer mean of the square of the 


vorticity 


” ~ “ ifj hemispheric moan of the terrain 


pressure 


F - i 2 n 

~ N ^i,j/ hemispheric mean of precipitation 


NOTE: N = ClCU-ICL+1) * (JMAX- 


2) 


rate 


)/2Ax]2, 


t)2Ax]2^ 
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Subroutine DIAGNL 


q* 


Subroutine DIAGNL computes and prints various 
2onal band diagnostic quantities in order to monitor the 
forecast. The zonal bands are ten degrees wide and range 
from the equator to the north pole. Table III-13 defines 
the diagnostic quantities computed by DIAGNL. 

There are four entry points in DIAGNL: IDIAGL, 
which initializes the diagnostic computations; DIAGLM, 
which computes all quantities except the eddy kinetic energy; 
DIAGLP, which computes the eddy kinetic energy; and DIAGPL, 
which prints the zonal band diagnostics at the end of the 
time step. 


r. Function POW 

. Y 

Functxon POW computes the power function, X , 

to six-digit accuracy to achieve a savings in computer time 

and is used to compute ir and in the heating routine. 

s . Subroutine SINCOS 

Subroutine SINCOS computes the sine and cosine 
of an argument to six-digit accuracy to achieve a savings 
in computer time and is used in the diagnostic routine DIAGNL. 


III-68 



‘4 


TABLE III-13 : DIAGNOSTIC QUAI'ITITIES COMPUTED BY 

SUBROUTINE DIAGML 


V ^ 4 i^cosX • ■ ■*4* V* • 1 sxriX * , m 


zonal mean of 


the tangential component of the wind 


^ -i ^sinX . V* 4 ,cosX. ,, zonal mean of 
"xt i,j ^t3 irJ/3^ 1/3 

the normal component of the wind 


mean of the kinetic energy 


zonal mean of the eddy kinetic energy 

Tt V “ ^ ^ 4 , , zonal mean of the temperature 

i/j 1/3/*^ 

^ ^ S |wj^ j y^, zonal mean of the vertical velocity 

L X f j ' 

“^L,K - ""i+l,j,k - “i-l,j,k + ''i.j+l.k - ■'i,j-l,k>''24«l^ 

zonal mean of the square of the divergence 

zonal mean of the square of the vorticity 


NOTE: L denotes the zonal band and N^ denotes the number 

Jb 

of grid points lying in that band 
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4. Output Overlay - Program OUTPE 

The function of the output overlay is to take the 
results of the forecast contained on file lOLD and to pro- 
duce meteorologically meaningful output; that iS/ to create 
surface pressure, precipitation and geopotential heights, 
temperatures and winds at standard levels for viewing by 
meteorologists. In addition, the output overlay writes a 
restart data tape on each output cycle (nominally set at 
twelve- hour intervals) . 

Figure III-13 shows the structure of the output over- 
lay, giving the names of all of the subprograms which 
constitute this overlay. The function of each of these 
routines is described in the sections that follow. 

Table III-14 defines the fields which are produced by 
the output overlay and written to TAPE4. 

Table III-15 shows the common blocks local to this 
overlay. 

a. Subroutine SQRTQUT 

Since the output of the forecast contained on file 
lOLD is data arranged in vertical slice format and the 
required output is data on horizontal planes, an interface 
routine to properly structure the data is required. Sub- 
routine SORTOUT takes forecast data in vertical slice format 
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TABLE III-14: CONTENTS OP THE OUTPUT DATA FILE - 

TAPE 4 PER OUTPUT CYCLE 

Record Data Field 

1 Forecast Surface Pressure/ Ps (mb) 

2 Precipitation (cm) 

3 Geopotential Heights of 1,000 mb Surface (meters) 


4 


ri 



11 

TI 

950 

It 

II 


II 

5 


ir 



ir 

It 

900 

It 

II 


II 

6 


ti 



IT 

|1 

850 

II 

II 


It 

7 


ti 



II 

It 

700 

II 

II 


II 

8 


It 



It 

11 

500 

M 

11 


II 

9 


ti 



IT 

IT 

400 

II 

11 


11 

10 


IT 



11 

ir 

300 

II 

II 


II 

11 


U 



11 

11 

250 

II 

IT 


11 

12 


ir 



II 

II 

200 

II 

II 


II 

13 


11 



It 

II 

150 

II 

11 


II 

14 


11 



IT 

11 

100 

TI 

11 


IT 

15 

Temperatures on 

the 1/000 

mb Surface 

( 

^C) 


16 


n 


ri 

11 

950 

II 

ri 


II 


17 


IT 


II 

II 

900 

II 

TI 


II 


18 


II 


rr 

ir 

850 

rt 

TI 


II 


19 


IT 


ri 

It 

700 

II 

11 


II 


20 


II 


II 

IT 

500 

II 

11 


tl 


21 


IT 


ti 

It 

400 

II 

11 


IT 


22 


II 


II 

II 

300 

It 

II 


II 


23 


II 


II 

It 

250 

IT 

IT 


It 


24 


II 


IT 

II 

200 

It 

II 


11 


25 


11 


II 

Tt 

150 

II 

II 


tt 


26 


IT 


It 

ir 

100 

II 

11 


IT 


27 

U 

winds 

on 

1000 

mb 

surface 

(m/sec) 




28 

V 

II 

IT 

1000 

It 

TI 


fl 




29 

U 

II 

IT 

950 

It 

II 


II 




30 

V 

IT 

IT 

950 

It 

II 


IT 
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TABLE III-14; OUTPUT DATA FILE - TAPE4 (Continued) 


Record Data Field 


31 

U 

winds 

on 

900 

mb 

surface 

(m/ sec) 

32 

V 

tt 

It 

900 

II 

IT 

II 

33 

U 

\t 

tt 

850 

II 

II 

Tl 

34 

V 

II 

tt 

850 

IT 

It 

II 

35 

U 

ir 

IT 

700 

II 

It 

IT 

36 

V 

II 

IT 

700 

tt 

It 

IT 

37 

U 

M 

tt 

500 

11 

11 

It 

38 

V 

tt 

11 

500 

II 

ft 

11 

39 

U 

II 

II 

400 

II 

II 

II 

40 

V 

Tt 

II 

400 

II 

11 

II 

41 

U 

II 

II 

300 

It 

II 

Tt 

42 

V 

II 

II 

300 

II 

II 

tl 

43 

U 

II 

11 

250 

II 

Tl 

It 

44 

V 

11 

II 

250 

11 

IT 

It 

45 

U 

tt 

II 

200 

ri 

ft 

ir 

46 

V 

ft 

If 

200 

it 

II 

tt 

47 

U 

11 

tl 

150 

IT 

11 

IT 

48 

V 

ir 


150 

II 

IT 

tr 

49 

U 

ir 

II 

100 

11 

IT 

II 

50 

V 

11 

ir 

100 

II 

It 

11 


NOTE: Each data field has an identification block of twenty 

words appended to it, which uniquely identifies the 
forecast base and forecast time plus parameters 
for plotting of the fields by the graphics package. 


i 

f 


1 ! TABLE III-15: 

f 


i ^ COMMON /WORK/ 

PP1(63,63,12) 

r ■ ; 

i i 

) i 

COMMON /BUFFP/ AO (3989) 


COMMON BLOCKS LOCAL TO 
THE OUTPUT OVERLAY - 
MODELS PECHCV AND PECHPV 


Used as temporary storage 
for the forecast fields on 
horizontal planes. 

Used as input/output buffer 
area for file lOLD. 


I 
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on sigma surfaces and constructs data on standard pressure 
surfaces and stores it on the random access file TAPE2. 

In order to construct a data field for output 
purposes, the forecast values must be interpolated or 
extrapolated to standard pressure levels as the file lOLD is 
swept, row by row. This process creates horizontal plane 
data fields. The interpolation process is straightforward 
using logarithmic interpolation in pressure to generate 
output data on presstire surfaces. However, when a pressure 
surface is below the terrain, simple extrapolation can give 
very unrealistic results and a more sophisticated proce- 
dure is required. 

The method described below is a variation of the one 
in use at FNWC and is based on the work of Bellamy (1945) . 

It basically uses the standard heights and temperatures of 
pressure surfaces and modifies them according to the condi- 
tions at a = 1.0. Define; 


1.0 


, 288.16 , ,, , TT , ^5.2561^ ^ 

^ ^ ^ T m o c / J 


6.5x10 


‘1013.25^ 


T 


1.0 


288.16-Z^ q(6.5x10 


- 1 


SFC 


P — ) 

, 288.16 , „ , s , ^5.2561^ t 

^ —3 oc/ J 


6.5x10 


'1013.25' 


[III-9] 
[III. 101 

[III. 11] 
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and 




[III. 121 


[111.13] 


Then if the difference 2^^ q - is less than 500 m. ;■ 

the height and temperature at standard level L is given by: 

ML) = ®l.o’ [III. 141 

and 

T(L) = Tg^od.) [1 + [III. 15] 

Otherwise, by 

Z(L) = tZsjn(L) - [1 + |(3S - Sj__(,)l [III. 161 

and 

T(L) = Tg^j^(L)[l + S] [III. 17] 

where Zg^jj{L) and Tg^p(L) are the standard height and tem- 
perature, respectively, of the standard pressure level, L. 

For graphical presentation, the surface pressure, geo- 
potential heights and temperatures are then smoothed to 
remove small-scale computational noise from the forecast 
values. The forecast fields are then written to the 
random access file TAPE2. Using subroutine PLTWRT which 
prefixes a twenty-word identification record to the fields, 
the forecast fields are written to the forecast file TAPE4. 
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b. 


Subroutine PLTWRT 


Subroutine PLTWRT prefixes a twenty-word identi- 
fication record to a data field and then writes the field 
to a specified output unit. The identification record 
contains forecast and fxeld identifiers and contouring 
information for the. graphical display package. 

c. Subroutine FILTR 

Subroutine FILTR is a non-linear smoothing routine 
used to remove computational noise from selected output 
forecast fields. Three types of filter can be chosen: a 

medium wave filter with a cutoff wave length of six grid 
lengths? a short wave filter with a cutoff wave length of 
four grid lengths; and a long wave filter with a cutoff 
wave length of thirty grid lengths. 

d. Subroutine SMOOTH 

Subroutine SMOOTH is a diffusive type of data 
field smoother that only smooths points whose latitude 
is less than thirty degrees. These points are then smoothed 
according to the formula 

P[ , = F. . + 0.1 (P) [III. 181 

1 1 J-/ 3 


Multiple passes are made over the data field. 



e. 


Subroutine PRT 


Subroutine PRT is a specialized printer contour 
routine which contours 63 x 63 data fields on four pointer 
pages. It is used to output a selected few data fields 
to get a quick look at the quality of the forecast. 

f . Subroutine RSTRTD 

Subroutine RSTRTD v;rites restart data to TAPE3 
during each output cycle for possible restart use in the 
event of computer malfunction or if it should be desired 
to continue a forecast at a later date. The format of 
TAPE3 is given in Table III-IO, which appears in Section 
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Code Structure - PEFHCV and PEFHPV 


This section describes the program structure of the 
forecast models PEFHCV and PEFHFV. It also defines the 
various files used by the programs, defines the common blocks 
that are used and describes each of the subprograms that 
constitute the models. 

The codes are written in CDC Fortran, which is a 
dialect of FORTRAN IV, for the CDC FORTRAN Extended Compiler 
C0PT=2 level) . The txvo programs, PEFHCV and PEFHFV, have 
been designed to use rotating mass storage along with central 
memory and either ECS or LCM. 

The forecast models are fairly large codes. In fact, 
the models consist of a set of two programs . PEPREP sorts 
out the analysis data into a vertical slice-pressure surface 
format. The main forecast model code consists of PEFHCV or 
PEFHFV. Section 1 describes the analysis data program, 
PEPREP, giving a block diagram of the program structure, 
and describing each of the subprograms. 

The forecast models are fairly large codes. Since 
the various sections of the codes — initialization, inte- 
gration and output — are each fairly large, the main 
forecast model codes have been overlaid to more efficiently 
use central memory. Section 2 describes the control overlay 
and gives a block diagram of the overlay structure, defines 
the I/O file structure and defines the common blocks that 
are global to all overlays . 




■ 
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Section 3 describes the initialization overlay defining 
the coiranon blocks local to this overlay and describing each 
of the subprograms that are used by this overlay. Similarly, 
Sections 4 and 5 describe the integration and output overlays, 
respectively. 

1. Analysis Data - Program PEPREP 

Due to the large size of the data base for the 187 x 187 
forecast models (approximately one million words) , it is 
much more efficient to sort out the analysis data into a 
pressure surface - vertical slice format prior to doing the 
sigma surface interpolation in the initialization section 
of the forecast model. Program PEPREP takes analyzed data 
on pressure surfaces and geophysical data, computes moistures 
and 50 mb winds, and, finally, sorts out this large mass of 
data producing a vertical slice - pressure surface format 
data file. 

The structure of the program PEPREP is shown in Figure 
III-14 and the various files used by PEPREP are defined in 
Table III-16. In addition, the program uses both rotating 
mass storage and either ECS or LCM (approximately 250,000 
words) as temporary data storage. The contents of the 
input data file, TAPEl, have been defined previously in 
Table III-8 in Section III-2 with the exception that the 
187 X 187 analysis data records are packed 2:1 prefixed by a 
20 word identification header. 



‘ ■ ■ ■ ■ ■■ - 
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TABLE III-16: 


FILE STRUCTURE OF ANALYSIS 
DATA PROGRAM - PEPREP 


OUTPUT 
TAPEl 
TAPE 2 

TAPES, 
TAPE 6 


Used to write various diagnostic messages 

Input data file generated by analysis program 

Random access file used for intermediate 
storage of data fields 

(JOLD, JNEW) Vertical slice analysis data 
on pressure surfaces for input to the main 
forecast models PEFHCV or PEFHFV 
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3 .* 


Subroutine TREAD 


Subroutine TREAD acts as a generalized input data 
routine to obtain analysis data from TAPEl. A storage 
array address and the data field ID are input to TREAD from 
the calling routine. TREAD then searches a table (array ITP) 
to find the location of the desired data field and reads the 
data field from TAPEl into the storage array. Various data 
field orderings on TAPEl may be accommodated simply by 
changing the information stored in array ITP. 

b. Subroutine DATAIO 

Subroutine DATAIO is a generalized routine for 
writing and reading data fields to and from the random 
access file TAPE2. Two entry points exist for this routine: 
DATAOUT, which uses the system routine WRITMS for writing 
to the filer and DATAIN, which uses the system routine READMS 
for reading from the file. 

c. Subroutine REPACK 

Subroutine REPACK repacks analysis data in a 2:1 
format for use in generating the data records for the vertical 
slice - pressure surface file. For example , if one has 
two packed analysis data arrays consisting of the u and v 
wind components r then REPACK puts these two arrays into one 
187 X 187 packed array with the u component of the wind in 
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the upper half of the word and the v component in the lower 
half of the word. The routine is also used for repacking 
other data fields. 

d. Subroutine MOIST 

Subroutine MOIST computes the relative hxaraidity 
on a pressure surface using the procedure based on the 
relative vorticity as described in Section II-E-6 . Input 
to the routine consists of the geopotential heights of the 
pressure surface. The routine then computes the relative 
vorticity on the pressure surface and, using the empirical 
relation given by Equation III. 133], computes the relative 
humidity on the pressure surface. Program PREPREP calls 
this subroutine five times to generate moistures on the 
1,000, 900, 700, 500 and 300 mb surfaces. 

e . Subroutine WIND 

Subroutine WIND computes the geostrophic wind 
components at 50 mb for use in the initialization of the 
winds. Input to the routine consists of the geopotential 
heights of the 50 mb pressure surface. 



2 . 


Main Overlay - Programs PEFHCV and PBFHFV 


The main overlay of the forecast model progams is 
very simple, since it only performs the function of declar- 
ing three global common blocks, defining various control 
variables and calling the actual working overlays. The 
overlay structure of the programs is shown in Figure III-15. 

The various files used by the programs are defined in 
Table III-17. Tables III-18 and III-19 define the variables 
that appear in the global common blocks CNTRL, and INDEX, 
respectively. Common block BUFPOINT consists of the offsets 
for the variables and contains the integer variables listed 
in Table III-3. 


a. Stibroutine BUFERR 

Subroutine BUFERR prints a diagnostic message 
when an error occurs on a BUFFER (IN or OUT) operation that 
tells where the error occurred. 
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FIGURE III-15: OVERLAY STRUCTURE OF FORECAST MODELS 

PEFHCV AND PEFHFV 








TABLE III-17: 


FILE STRUCTURE OF 
FORECAST MODELS 
PEFHCV AND PEFHPV 


OUTPUT 
TAPE 2 
TAPES 

TAPE 4 

TAPES 

TAPES 

TAPES 


Used to write various diagnostic 
messages 

Random access file used by output 
overlay (four records only) 

Restart data tape; contains a 
header record and contents of 
current lOLD file 

Output data file; contains fore- 
cast fields to be selectively 
plotted by graphics program and 
for use by the analysis programs 

Input data file for initialization 
overlay; generated by program PEPREP 

(lOLD, INEW) integration data 

(lOLD, INEW) integration data 


TABLE III-18: DEFINITION OF VARIABLES IN COMMON 

BLOCK /CNTRL/-MODELS PEPHCV AND PEFHFV 


lOLD 

INEW 

ISW 


ID (2) 

LIMI 

LIM 

ITAU 

ITSPHR 

NN 


NNOUT 

LCO 

L12 


Pile name of integration data (TAPE8 or TAPE9) 

Pile name of integration data (TAPES or TAPE 9) 

=1, SSWl is "on", signifying normal start xi^ith 
data from TAPEl 

~ 2 , SSWl is "off", signifying restart with data 
from TAPES 

1st word, forecast identifier (date-time group) 

2nd word, forecast data field identifier 

Starting cycle 

Ntimber of integration and output cycles 
(2 X number of days + 1) 

Forecast time indicator (hours) 

= 0 during initialization and first output cycle 
= 12*LCO for normal integration cycle 

Number of integration time steps per hour 

Time step counter, reset after each output cycle, 
also NN=0 signifies a Matsuno time step 

Number of time steps between normal output cycles, 
usually equal to 12 x ITSPHR 

Output cycle counter, updated by 1 each output cycle 

Indicator set at 12-hour output cycle, used to 
zero rainfall 


DAY 

HR 

DT 

TDT 

MATS 


MFRQ 


Number of days from January 1 of the forecast 

Hour of the day, GMT 

Time step. At 

Twice the time step, 2 At 

Logical variable used during Matsuno time step 
= .FALSE, during first half 
= .TRUE, during second half 

Number of hours between Matsuno time steps 
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TABLE 

ICL 

ICU 

IMAXM 

IMAX 

JMAXM 

JMAX 

KMAXQ 

KMAXM 

KMAX 

LINE 

NPLD 

NFLDEX 

NVARS 

NVARL 


III-19: DEFINITION OF VARIABLES IN COMMON 

BLOCK /INDEX/ -MODELS PEFHCV AND PEFHFV 

Start of computation on data row (=2) 

End of computation on data row (=IMAX-1) 

Length of data row -1 

Length of data row 

Number of data rows -1 

Number of data rows 

Number of moisture levels 

Number of vertical levels -1 

Number of vertical levels 

Number of row being worked on 

Size of a data record on TAPE2 (63 x 63) 

Expanded data record (187 x 187) 

Size of the AJM and AJP buffers 

Size of the AJ buffer and a data record on 
TAPE 8 or TAPE9 


3 


Initialization Overlay - Program INITEX 


The function of the initialization overlay is to taj^e 
as input data the analysis data in vertical slice - pressure 
surface format generated by Program PEPREP and to generate 
an initial state for the integration overlay which actually 
performs the forecast. That is, the interpolation of data 
fields to the PEM sigma surfaces as described in Section II-E 
is performed by the initialization overlay. 

Figure III-16 shows the structure of the initialization 
overlay, giving the names of the subprograms which consti- 
tute the overlay. 

Program INITEX performs the simple executive function 
of either calling overlay RSTART to initialize from a pre- 
viously written restart file, or by calling overlay SORT to 
initialize from the analysis data file written by Program 
PEPREP . 


a. Program SORT 

Program SORT takes the analysis data in vertical 
slice - pressure surface format as generated by program 
PEPREP and generates the integration data file in vertical 
slice - sigma surface format. The interpolation to sigma 
surfaces proceeds one vertical slice at a time sweeping the 
input data file and creating the integration data file. The 
interpolation process for converting pressure surface data 
to sigma surfaces is described in Section II-E. 
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FIGURE I II- 16 : INITIALIZATION OVERLAY STRUCTURE 

MODELS PEFHCV AND PEFHFV 




b. Program RSTART 


Program RSTART reads the restart data file/ TAPES/ 
when a forecast run is being restarted and initializes the 
forecast model in order to continue the forecast. Restarts 
are controlled by sense switch 1, If SSWl is "ON”/ then 
a normal run is being made; if SSWl is "OFF"/ then a restart 
run is being made. 

A restart data tape is written during each output 
cycle so if the computer goes down/ or if for some reason/ 
a forecast is to be continued at a later time/ the forecast 
can be resumed from the point at which the last restart tape 
was written. Table III-IO gives the format of the restart 
data tape TAPES . 

TnuB/ program RSTART reads in the header record 
from TAPES and initializes the required control variables. 
Next the routine reads in the forecast data records and 
recreates file lOLD. The forecast rim will now proceed 
from this point in the forecast. 
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4. 


Integration Overlay - Program INTGPE 


The function of the integration overlay is to numer- 
ically integrate the primitive equations; that is, to step 
the data fields forward in time to produce a forecast. 

Input to this overlay consists initially of the data field 
specification produced by the initialization overlay and, 
subsequently, the current state of the forecast contained 
on file lOLD. After initialization, a forecast consists of 
cycling through the integration and output overlays until 
the desired forecast time has been reached. 

Figure III-IO in Section III-B-3 shows the structure of 
the integration overlay, giving the names of all the subpro- 
grams which constitute this overlay. Table III-20 defines 
the variables that appear in the common blocks that are 
local to the integration overlay for programs PEFHCV and 
PEFHFV. Table III-ll in Section III-B-3 defines the varia- 
bles in common block /CON/ which is the same for all four 
forecast model programs. 

With the exception of the indexing method used to access 
variables in the working arrays, the subprograms in the inte- 
gration overlay are identical to those occurring in models 
PECHCV and PECHFV. Therefore, the descriptions of the sub- 
programs are not repeated in this section and the reader is 
referred to Section III-B-3. 


TABLE III- 20; DEFINITION OF VARIABLES IN 

COMMON BLOCKS LOCAL TO THE 
PROGRAMS PEFHCV AND PEPHFV 


Common Block /LCM/ 
AIN(NVARL,3) 

Common Block /BUFFS/ 
AJM(NVARS) 

AJ (NVARL) 

AJP (NVARS) 

Common Block /LOCAL/ 
PIKAP(187) 
W(187,KMAXM) 

DUM(187,KMAX+4) 

a) PILN{187,3) 

b) UPMC187) 

VPM(187,3) 

AB(187,KMAX) 

c) PIP (187,3) 

PEPS (187) 


Circular input buffer for rows JM, 

J and JP; resides in either ECS or LCM 


Buffer for partial row JM in central 
memory 

Buffer for full row J in central memory 

Buffer for partial row JP in central 
memory 


on row J 

Vertical velocities for all levels on 
row J 

Dummy array 

NOTE: DUM is equivalenced to the 

following variables 

Jin (tt) on rows JM, J and JP 

on row J at a particular level 

V7T 

on rows JM, J and JP at a particular 

level 

2 

Aqm (5a+6B) at all levels on row J 
2Ax 

TT reduction factor to sea level 

ir change due to nonlinear smoothing 
operator 
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TABLE III-20: Definition of Variables (Continued) 


Common Block /LTEND/ 
TENDP(187,2*KMAX) 
DPIDX(187,KMAX) 

DPIDY (187,KMAX) 

TCOR(187,3,mAX) 


a) PIT(1S7) 
TT(187,KMAX) 

QT(187,KMAXQ) 

b) UT(187,KMAX) 
VT (187,KMAX) 

PC(187,KMAXQ) 


Dummy tendency array 

Pressure-gradient— force terms in x 
direction for all levels on roir? J 

Pressure-gradient— force terms in y 
direction for all levels on ro\^ J 

Temperatures corrected for vapor 
pressure at all levels on rows JM, 

JP 

/LTEND/ resides in central memory 
on the Cyber 175 and in LCM on the 
Cyber 76 

TENDP is equivalenced to the 
following variables: 

terrain pressure tendency on row J 
tendencies in irT for all levels on 


tendencies in irg for moisture levels 
on row J 


J and 
NOTE ; 

NOTE : 

3t ^ 

3TrT 
3t " 
row J 

3t " 


— ^ , tendencies in ttu for all levels on 
row J 

37T V * 

, tendencies in irv for all levels on 
row J 

NOTE ; TCOH is equivalenced to the 
following variables: 


Pressure at moisture levels on row J 
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Definition of Variables (Continued) 


TABLE III-20; 


QC(187,KMAXQ) 

QCS(187,KMAXQ) 

TC (187,KMAX) 
RP(187) 


Vapor pressure at moisture levels on 
row J 

Saturation vapor pressure at moisture 
levels on row J 

Temperature at all levels on row J 
Rainfall on row J 


NOTE: The following common block only appears in the 

model PBFHCV. 

Common Block /SCM/ 

IBDF (NVARL) Input buffer for file lOLD 

OBUP(NVARL) Output buffer for file INEW 


NOTE: Symbolic dimensions are given where appropriate since 

the size of the array varies with the forecast model. 
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5. 


Output Overlay - Program OUTPE 


The function of the output overlay is to take the 
results of the forecast contained on file lOLD and to 
produce meteorologically meaningful output; that is, to 
create surface pressure, precipitation and geopotential 
heights , temperatures and winds at standard levels for 
viewing by meteorologists. In addition, the output over- 
lay writes a restart data tape on each output cycle (nomi- 
nally set at twelve-hour intervals) . 

Figure III-17 shows the structure of the output overlay, 
giving the names of all of the subprograms which constitute 
this overlay. The function of each of these routines is 
described in the sections that follow. 

Table III-14 in Section III-B-4 defines the fields which 
are produced by the output overlay and written to TAPE4. 

Table III-21 shows the common blocks local to this overlay. 

a. Subroutine SORTOOT 

Since the output of the forecast contained on 
file lOLD is data arranged in vertical slice format and 
the required output is data on horizontal planes, an inter- 
face routine to properly structure the data is required. 
Subroutine SORTOUT takes forecast data in vertical slice 
format on sigma surfaces and constructs data on standard 


pressure surfaces. 


















TABLE III-21: COMMON BLOCKS LOCAL TO THE OUTPUT 

OVERLAY - MODELS PEFHCV AND PEPHPV 


COMMON /WORKX/ DUMB (20), 
FA(X87,X87) 

COMMON /BUPFP/ AO (3969) 
COMMON /LCM/ PL (X87 ,X87, 3) 


Used as temporary storage 
for the forecast fieXds on 
horizontaX pXanes and as an 
input buffer area for fiXes 
lOLD and the anaXysis data 
fiXe. 

Used as a working storage 
area. 

Buffer area in ECS or LCM 
used during the interpoXation/ 
sorting process. 


In order to construct a data field for output 
purposes, the forecast values must be interpolated or 
extrapolated as explained in Section to standard 

pressure levels as the file lOLD is swept, row by row. This 
process creates horizontal plane data fields. For graphical 
presentation, the surface pressure, geopotential heights and 
temperatures are then smoothed to remove small-scale compu- 
tational noise from the forecast values. Using subroutine 
PLTWRT which appends a twenty-word identification record to 
the fields, the forecast fields are x-i^ritten to the forecast 
file, TAPE4. 

b. Subroutine PLTWRT 

Subroutine PLTWRT prefixes a twenty-word identi- 
fication block to a data field, packs the data field into a 
2;1 format, and then writes the packed plot field to a 
specified output unit. The identification block contains 
forecast and field identifiers and contouring information 
for the graphical display package, 

c. Subroutine FILTR 

Subroutine FILTR is a non-linear smoothing routine 
used to remove computational noise from selected output 
forecast fields. Three types of filter can be chosen; a 
medium wave filter with a cutoff wave length of six grid 
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lengths; a short wave filter with a cutoff wave length of 
four grid lengths; and a long wave filter with a cutoff 
wave length of thirty grid lengths. 

d. Subroutine SMOOTH 

Subroutine SMOOTH is a diffusive type of data 
field smoother that only smooths points whose latitude 
is less than thirty degrees. These points are then smoothed 
according to the formula 

FI . = F. . + 0.1 (F) [III. 19] 

1/ J Ir J 

Multiple passes are made over the data field. 

e. Svtbroutine PRT 

Subroutine PRT is a specialized printer contour 
routine which contours 63 x 63 data fields on four printer 
pages. It is used to output a selected few data fields to 
get a quick look at the quality of the forecast. 

f . Subroutine RSTRTD 

Subroutine RSTRTD \«rites restart data to TAPE3 
during each output cycle for possible restart use in the 
event of computer malfunction or if it should be desired 
to continue a forecast at a later date. The format of 
TAPE3 is given in Table III-lO, which appears in Section 
III-B-2. 
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Subroutine AVE 


Subroutine AVE transforms forecast data from a 
187 X 187 grid to a 63 x 63 grid so that contours may be 
plotted using subroutine PRT. Every third point along the 
boundaries of the 187 x 187 grid is moved to the boundaries 
of the 63 X 63 grid. The interior points are averaged via 
the following algorithm; 


16 


( 4F. . + . + F. + F. , . + F. . t) 

1,: 1+1,3 1,3+1 1-1,3 1,3-1 


^i+l,j-l ^i+l,j+l ^i-l,j+l '*■ ^i-l,j-l ) [III. 20] 
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